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1  Introduction 


Breast  cancer  is  the  most  frequently  diagnosed  malignancy  among  women  in  the  United 
States.  In  1999,  the  American  Cancer  Society  estimated  that  175,000  women  would  be 
newly  diagnosed  with  breast  cancer  and  that  43,300  would  die  from  the  disease  [1].  Breast 
cancer  accounts  for  29%  of  all  cancers  detected  and  16%  of  all  cancer  deaths,  and  ranks  as 
the  second  leading  cause  of  death  from  cancer  among  women  in  the  United  States  [1].  Five 
year  survival  rates  are  generally  very  high  (93%)  for  breast  cancer  staged  as  being 
localized,  falling  to  72%  for  regional  disease  and  only  18%  for  distant  disease  [2],  The  early 
detection  of  breast  cancer  is  clearly  a  key  ingredient  of  any  strategy  designed  to  reduce 
breast  cancer  mortality. 

The  goal  of  this  project  was  to  develop  computerized  tools  that  will  refine  the  perception  of 
mammographic  features  (including  lesions,  masses  and  calcifications).  Our  research  efforts 
were  geared  towards  improving  the  local  mammographic  viewing  environment  by 
selectively  processing  mammograms  for  presence  of  different  features,  and  towards 
providing  a  better  global  mammographic  viewing  environment  by  fusing  together  locally 
processed  sections  of  images.  By  improving  the  visualization  of  breast  pathology,  the 
chances  of  early  detection  of  breast  cancers  can  be  increased  (quality  improved)  while  less 
time  to  evaluate  mammograms  for  most  patients  required  (costs  lowered). 

We  were  investigating  a  methodology  for  accomplishing  mammographic  feature  analysis 
through  multiscale  representations.  In  this  report,  we  present  a  scheme  for  local 
enhancement  and  fusion  of  clinically  significant  features.  We  devised  a  wavelet  transform 
that  is  flexible  enough  for  incorporation  of  a  variety  of  enhancement  methods  and  used  the 
derived  wavelet  framework  for  enhancement  of  microcalcifications,  circumscribed  masses, 
and  stellate  lesions. 

In  the  following  sections,  we  briefly  overview  the  contents  of  the  report,  list  publications, 
and  explain  notation  that  we  use. 

1.1  Overview  of  Contents 

Wavelet  transform  forms  the  framework  of  our  contrast  enhancement  technique,  and 
Section  2.1  is  dedicated  to  it.  First,  Section  2.1.1  motivates  the  use  of  redundant 
representations  instead  of  orthogonal  and  biorthogonal  wavelet  transforms.  The  transform 
itself  is  then  described  in  Section  2.1.2,  while  Section  2.1.3  deals  with  fast  transform 
implementations  issues. 

Mammographic  image  enhancement  methods  are  typically  aimed  at  either  improvement  of 
the  overall  visibility  of  features  or  enhancement  of  a  specific  sign  of  malignancy.  Section  2.2 
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presents  a  synthesis  of  the  two  paradigms  by  means  of  image  fusion.  In  Section  2.2.1,  the 
suitability  of  steerable  dyadic  wavelet  transform  for  image  fusion  was  compared  to  two 
transforms  popular  in  fusion  applications,  the  gradient  pyramid  and 
orthogonal/biorthogonal  wavelet  transform,  by  means  of  both  subjective  and  objective 
criteria.  Section  2.2.2  next  reports  on  local  enhancement  strategies  targeting 
microcalcifications,  circumscribed  masses,  and  stellate  lesions. 

Evaluation  of  our  method  is  presented  in  Section  2.3.  The  method  was  compared  to 
histogram  equalization  and  unsharp  masking  on  image  phantoms  in  Section  2.3.1  and  on 
phantoms  blended  into  mammograms  in  Section  2.3.2.  The  comparisons  were  performed 
using  objective  measures,  but  behavior  of  the  three  techniques  was  assessed  from  the 
perceptual  point  of  view  as  well. 

1.2  Notation 

We  use  symbols  N,  Z,  and  R  for  the  sets  of  naturals,  integers,  and  reals,  respectively. 
L2(R)  and  L2(R2)  denote  the  Hilbert  spaces  of  measurable,  square-integrable  functions 
f[x )  and  f(x,y),  respectively. 

The  inner  product  of  two  functions  f(x)  G  L2{R)  and  g(x)  G  L2{R )  is  given  by 

/OO 

fix)  gix)  dx. 

-OO 

The  norm  of  a  function  fix )  G  L2iR )  is  defined  as 

ll/ll  = 

The  convolution  of  functions  fix )  G  L2(i? )  and  gix)  G  T2(K)  is  computed  as 

/OO 

fit)  gix  -  t)  dt, 

-OO 

and  the  convolution  of  two  functions  fix,y)  G  L2(i?2)  and  gix,y)  G  L2(jR2)  equals 

/OO  POO 

/  fitx,ty)gix-tx,y-ty)dtxdty. 

-oo  J  — OO 

The  Fourier  transform  of  a  function  fix)  G  L2iR)  is  defined  as 

~  f  oo 

/(w)  =  /  /(a;)e~JWI  dx, 

J  —  OO 

and  the  Fourier  transform  of  a  function  fix,y)  G  L2(Jt 2)  is  equal  to 

fiux,LUy)=[  [  fix,y)e~j{-WxX+UyV')dxdy. 

J  —  oo  J  — OO 
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l2(Z)  and  l2(Z2)  stand  for  the  spaces  of  square-summable  discrete  signals  f(n)  and 
f(nx,ny),  respectively. 

The  ^-transform  of  a  discrete  signal  f(n )  G  12{Z)  is  defined  as 

OO 

F(z)=  Y.  /(»)2_“- 

n=— oo 

The  convolution  of  discrete  signals  f(n)  G  l2 (Z)  and  g(n)  G  12{Z)  is  equal  to 

OO 

f*g(n)=  Y  f(m)g(n-m ), 

m=— oo 

and  the  convolution  of  discrete  signals  f(nx,ny)  G  l2(Z2)  and  g(nx,ny )  G  l2(Z2)  is  given  by 

00  00 

f*g{nx,ny)  =  Y  £  f(mx,my)g(nx-mx,ny-my). 

mx=— oo  my  — — oo 

The  Fourier  transform  of  a  discrete  signal  f(n)  G  l2(Z)  is  equal  to  the  ^-transform 
evaluated  on  the  unit  circle 

OO 

F{a>)=  E 

n=— oo 

and  the  Fourier  transform  of  a  discrete  signal  f(nx,ny)  G  12{Z2)  is  defined  as 

00  00 

F(ujx,ujy)  =  X)  E  /(«.,«„)  6-^"-^"-). 

tlx  = — OO  7ly  ——OC 

For  later  use,  we  define  the  following  functions: 


1.  the  unit  impulse  function 


2.  the  unit  step  function 


f  1  for  x  =  0 
1  0  otherwise, 

f  1  for  x  >  0 
|  0  for  x  <  0, 


3.  the  rectangular  function 


rect(x) 


1  for  |a;|  <  \ 
0  for  |a;j  >  |, 


4.  the  sine  function 


5.  the  unit  impulse  sequence 


sinc(a;) 


sin(7ra;) 

7TX 


and 


5(n)  := 


1  for  n  =  0 
0  otherwise, 


where  x  G  R  and  n  G  Z. 
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2  Body 

2.1  Wavelet  Transform 

Wavelet  transform  provides  the  framework  for  both  our  enhancement  and  fusion 
algorithms,  and  it  is,  therefore,  important  that  it  does  not  introduce  artifacts,  enables 
directional  multiscale  analysis,  and  can  be  implemented  efficiently.  In  Section  2.1.1,  we 
mention  problems  that  stem  from  the  lack  of  translation  and  rotation  invariance  of 
orthogonal  and  biorthogonal  wavelet  transforms.  Next,  in  Section  2.1.2,  we  describe 
approximations  of  steerable  wavelets  and  employ  the  multiscale  spline  derivatives  with 
both  first  and  second  derivative  wavelet  decomposition.  Our  annual  reports  contain  details 
on  higher  order  transform  derivation  and  implementation  with  reconstruction  ( i.e .,  inverse 
transform)  being  quite  cumbersome.  Here  presented  transform  allows  for  both  directional 
and  isotropic  processing  of  mammograms  while  being  shift-invariant  and  aliasing-free. 
Efficient  implementation  of  the  transform  is  highly  desirable,  and  Section  2.1.3  presents  a 
filter  bank  with  filter  implementations  taking  advantage  of  symmetry  and  antisymmetry  for 
faster  processing. 

2.1.1  Shortcomings  of  Traditional  Methods  of  Wavelet  Analysis 

Analyzing  images  across  multiple  scales  and  resolution  has  become  a  powerful  tool  for 
solving  compelling  problems  in  computational  vision,  image  processing,  and  pattern 
recognition.  Wavelet  theory  encompasses  multiscale  and  multiresolution  representations, 
such  as  subband  filtering  [3],  image  pyramids  [4],  and  scale  space  filtering  [5],  into  a  unified 
mathematical  framework.  In  the  area  of  image  processing,  there  remain  few  research  areas 
to  which  wavelet  analysis  has  not  been  applied.  For  example,  problems  in  image 
compression,  denoising,  restoration,  enhancement,  registration,  fusion,  segmentation,  and 
analysis,  have  all  been  approached  with  distinct  kinds  of  wavelet  processing. 

Though  ubiquitous,  wavelet  analysis  is  not  without  problems  of  its  own.  Lack  of 
translation  invariance,  one  of  the  major  problems  of  the  wavelet  transform  [6],  is  in 
multiple  dimensions  accompanied  with  lack  of  rotation  invariance. 

Wavelet  transform  in  its  most  commonly  used  orthogonal  or  biorthogonal  forms  is  not 
translation  and  rotation-invariant.  By  translation-invariant  transform,  we  mean  a 
transform  that  commutes  with  a  translation  operator.  Since  we  will  deal  primarily  with 
discrete  transforms  in  this  work,  we  constrain  the  translation  parameter  to  integer 
multiples  of  a  sampling  period. 

Lack  of  translation  invariance  of  the  discrete  wavelet  transform  is  illustrated  in  Figure  1. 
Here,  we  can  clearly  see  how  a  translation  of  the  input  signal  by  one  sample  results  in  a 
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completely  different  set  of  transform  coefficients  (orthogonal  wavelet  DAUB41  [6]  was  used 
in  this  experiment). 

Noninvariance  under  translations  of  an  orthogonal  and  biorthogonal  wavelet  transform  is 
due  to  lower  sampling  density  at  coarser  scales.2  A  straightforward  way  of  dealing  with  this 
problem  is  to  construct  a  redundant  transform  by  using  the  same  sampling  frequency  for 
the  input  signal  and  all  scales  of  the  transform.  A  filter  bank  implementation  of  such  a 
transform,  called  “algorithme  a  trous”  [7],  is  based  upon  the  fact  that  downsampling 
followed  by  filtering  is  equivalent  to  filtering  with  the  upsampled  filter  before  the 
downsampling,  as  shown  in  Figure  2. 

Lack  of  rotation  invariance  is  another  shortcoming  of  traditional  (i.e.,  orthogonal  and 
biorthogonal)  wavelet  techniques.  In  defining  rotation  invariance,  we  are  a  bit  less  strict 
than  with  translation  invariance.  We  do  not  require  that  the  transform  commutes  with  a 
rotation  operator  here.  Even  in  the  case  of  a  simple  filtering,  this  would  limit  us  to 
circularly  symmetric  filters  only.  Our  requirement  for  analysis  is  a  transform  that  enables 
rotation-invariant  processing.  As  an  example  of  such  a  transform,  let  us  consider  filtering 
with  the  first  derivative  of  a  two-dimensional  Gaussian  probability  density  function  in  two 
directions,  specifically,  along  x  and  along  y-axis.  By  linearly  combining  the  results  of 
filtering  in  these  two  directions,  filtering  with  the  first  derivative  of  a  Gaussian  in  any 
direction  can  be  computed.  This  fact  was  used  by  Canny  [8]  for  edge  detection.  A 
determined  edge  direction  rotates  as  an  input  image  is  rotated. 

After  choosing  the  fundamental  properties  of  the  transform,  one  must  decide  upon  the 
basis  functions  to  be  applied.  For  our  studies,  we  selected  basis  functions  that  well 
approximated  derivatives  of  a  Gaussian,  because  (1)  the  Gaussian  probability  density 
function  is  optimally  concentrated  in  both  time  and  frequency  domain,  and  thus  suitable 
for  time-frequency  analysis,  (2)  higher  order  derivatives  of  a  Gaussian  can  be,  similar  to 
the  first  derivative,  used  for  rotation-invariant  processing  [9],  and  (3)  the  Gaussian  function 
generates  a  causal  (in  a  sense  that  a  coarse  scale  depends  exclusively  on  the  previous  finer 
scale)  scale  space  [10].  The  last  property  makes  possible  scale-space  “tracking”  of  emergent 
features. 

xThe  number  in  DAUB4  refers  to  twice  the  order  of  the  wavelet  (i.e.,  two  in  this  case). 

2In  practice,  since  analysis  is  performed  over  a  finite  range  of  scales,  a  discrete  wavelet  transform  is 
translation-invariant  by  translations  determined  by  the  coarsest  scale  (e.g.,  sixteen  samples  for  the  analysis 
from  Figure  1)  [6]. 
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Figure  2:  Filter  bank  implementation  for  (a)  a  discrete  wavelet  transform  and  (b)  “algorithme 
a  trous”  decompositions  for  three  levels  of  analysis. 

2.1.2  Multiscale  Spline  Derivatives 

We  define  a  steerable  dyadic  wavelet  transform  of  a  function  s(x,y)  G  L2(R 2)  at  a  scale 
2m,  m  G  Z,  as  [11] 

W2ms(a;,y)  =  s*  ipl2m(x,y),  (1) 

where  tp2m(x,y)  denotes  ij)2m{x,y)  rotated  by  9{ ,  tp2m(x,y)  =  2~2mip(2~mx,2'~my),  ip{x,y)  is 
a  steerable  wavelet  that  can  be  steered  with  I  basis  functions,  and  6{  =  with 
i  €  {1, 2, ... ,  /}.  (For  introduction  to  steerability,  please  refer  to  our  first  annual  report.) 
Analogous  to  the  one-dimensional  case,  we  require  the  two-dimensional  Fourier  plane  to  be 
covered  by  the  dyadic  dilations  of  xj)l{2mux,2Tnu)y):  there  must  exist  A3  >  0  and  B3  <  oo 
such  that 

00  I 

A,  <  £  2<B3  (2) 

m=— oo  i=l 

is  satisfied  almost  everywhere. 

If  (nonunique)  reconstructing  functions  x1?™  (x,  y)  are  chosen  such  that  their  Fourier 
transforms  satisfy 

oo  I 

E  =  (s) 

m=— oo  i=l 

the  function  s(x,y)  may  be  reconstructed  from  its  steerable  dyadic  wavelet  transform  by 

oo  I 

s{x,y)=  *X2m(x,y),  (4) 

m=— oo  i— 1 

where  X2m{x,y)  denotes  X2m{x,  y)  rotated  by  9i  and  X2m(x,y)  =  2~2mx(^mx,  2~my). 
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We  choose  wavelets  that  are  steerable  analogs  to  the  one-dimensional  derivatives  of  central 
B-spline  wavelets  [12]: 

/si  n(^)\p+d+1 

1>{ur,w»)  =  (jurcos(ujo))d  y— ,  (5) 

where  wr  =  y/u>%  +  ujg  =  axg[pjx,uy),  and  d  G  {1, 2}.  These  wavelets  can  be  steered  with 
d+1  basis  functions. 

Wavelets  (5)  are  equal  to  d-th  order  derivatives  of  circularly  symmetric  spline  functions  in 
the  direction  of  x-axis  (note  that  knots  for  these  splines  are  circles).  To  implement  the 
transform  efficiently,  we  approximate  the  wavelets  with  x-y  separable  wavelets 

il>{x,y)  =  -~-^T~Pp+d(y)>  (6) 

where  /3p(x)  denotes  the  central  B-spline  of  order  p. 

Based  on  the  fact  that  B-splines  tend  to  a  Gaussian  probability  density  function  as  their 
order  increases,  it  is  easy  to  see  that  both  wavelets  (5)  and  (6)  converge  to  the  same 
functions  (i.e.,  d-th  order  derivatives  of  the  normalized  Gaussian  in  the  direction  of  x-axis) 
as  p  — >  oo.  In  order  to  steer  wavelets  t/>(x,  y)  given  by  (6)  (note  that  steering  will  be  only 
approximate,  since  these  wavelets  are  not  steerable),  we  need  to  find  basis  functions  that 
will  approximately  steer  ijj(x,y).  To  accomplish  this,  we  take  advantage  of  the  property  of 
circularly  symmetric  functions  that  rotations  of  their  directional  derivatives  are  equal  to 
directional  derivatives  in  rotated  directions: 

\ddQc{x,y)\  ddgc(x,y) 

00  I  dnd  }  dnd0  ’ 

where  7 Zg0  stands  for  rotation  by  60,  decg^  =  dVgc(x,y),  Qc(x,y)  is  a  circularly  symmetric 
function,  and  nSo  denotes  vector  ft  =  (cos  6,  sin  9)  rotated  by  9q. 

Let  us  choose 

g(x,y)  =  pp+d(x)Pp+d(y), 

which  is  approximately  circularly  symmetric  function  for  higher  order  splines.  A  rotation  of 
■0(x,y)  =  a  jjj-- l’ ^  from  Equation  (6)  by  90  can  therefore  be  approximated  by 


ip6o(x,y) 


ddg(x,y) 


=  E 


nd~lnl 


i  dd  iPP+d{x)  dl/3p 


dxd~l 


where  n  =  (cos  do,  sin  90)  =  ( nx,ny ).  (Note  that  in  case  of  Gaussian,  which  is  both  x-y 
separable  and  circularly  symmetric,  Equation  (7)  becomes  exact.) 
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To  derive  an  algorithm  for  the  fast  computation  of  the  transform,  we  introduce  two 
smoothing  functions  such  that 

oo  I 

fax,  Wy)  tp(u>x,  U)y)  =  Y,  £  2 mU)y)  Xl(2mWxi  2muy).  (8) 

m= 0  i=l 

Using  a  set  of  basis  functions  (7)  that  approximately  steer  wavelets  (6),  we  want  to 
construct  a  transform  such  that  Equations  (1)  through  (4)and  (8)  will  be  valid  (superscript 
i  must  be  viewed  now  as  an  index,  rather  than  rotation  by  0j). 

Let  F(uj )  be  a  digital  filter  frequency  response  and  let  us  denote 

=  e^sF(u) 

with  s  being  a  filter  dependent  sampling  shift  needed  to  obtain  finite  impulse  response 
(FIR)  filters. 

In  frequency  domain,  we  can  express  basis  functions  from  (7)  as 

4)i+1(u}x,ujy)  =  GdSsl(oJx)Gls(ojy)Pp+i(ujx)pp+d-i(ujy),  i  €  {0, 1, 2},  (9) 

where  Gd(u>)  is  given  by 

GV)  =  ^"*(2? sin  (!))',  (10) 

where  d  is  the  order  of  the  derivative,  d  G  {1,2},  the  sampling  shift  for  filter  (10)  is 
s  =  and  G°(u)  =  1. 

Since  we  are  interested  in  the  first  and  second  derivative  wavelets,  we  impose 


and  choose 

X  Fs[u)x)T[Uy)Pp(u)x)Pp-2{^y)i 

(11) 

X  (wx,wy)  Fs(u)x)Ks(u)y')0p_i(uix)j3p_i(u)y), 

(12) 

X  (^nwy)  T(yLOx)Ks(^UJy^pp_2{^x)Pp{^Jy)j 

(13) 

where 

*■'->  ■  iw  (?»“' (£  Mi))T'  ■ 

(14) 

with  the  sampling  shift  for  Kd(u )  being  the  same  as  the  one  for  Gd(oj),  and 

T(u,)  =  |ff(a>)|2 

(15) 
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with 


(16) 


H(u)  =  e>"  (cos  (!))  , 


the  sampling  shift  for  H(co)  being  s=  0>flbnod2, 

Using  the  relation  derived  in  the  previous  report 

4,(2u;)  =  H_s(u)/3P(  u)  (17) 

together  with  Equations  (9)  and  (11)  through  (13)  with  Equation  (8)  results  in 

G2{ujx)K2(u}x)T(ojy)  +G1(ojx)K1(ojx)Gl(uy)K1(uy)  +  T(ujx)G2(uy)K2(uy)+ 

+\H{u,x)H(uy)\2  =  1. 

Next,  we  derive  a  filter  bank  implementation  of  the  transform.  Assuming  a  bandlimited 
input  signal  s(ux,uy)  =  0  for  \ux\  >  tt  or  \u>y\  >  n  and  using  Shannon’s  sampling  theorem  in 
two  dimensions  [13]  with  Equation  (1)  and  basis  functions  from  Equation  (9),  we  can  write 


roc  rOO 

1  s(x,  y)=  /  s(4,  iy)  sine (tx  -  ix)  sine (ty  -  *„)• 

J —oo  J  — oo  -  ■  _  „ 


-oo  00  ix~~~oo  iy=—c>o 


9-s  imx)Pp+i{x  tx  wix)  9—s{my)Pp+d-i(y  ty  my)  dtx  dty, 

YYlx —  OO  TTly  — — "  OO 

where  i  G  {0, 1, 2}  as  in  Equation  (9). 

We  approximate  sine  functions  with  r-order  cardinal  splines,  then  use  the  relation  between 
cardinal  and  B-splines 

OO 

Vr(x)=  b~1(i)^r(x  —  i) 

i=  —  OO 

and  get 

W^2ms(x, y)  _  _  }  —  S(ujx,  Uy')  Br  (ujx)  Br  (k,j/)  •5p+r+i-i-i(k’i)’ 

x — nx  ,2/ — fty 

m—  1 

'Bp+r+d- i+1  (idy  )  GtJi 2mux)  Gls(2mujy)  [J  Hp+si(2ncjx)Hp_+sd~i(2nuy),  (18) 

n=0 

where  B~l{ui)  denotes  the  Fourier  transform  of  the  direct  B-spline  filter  of  order  p.  Table  1 
shows  the  ^-transforms  of  direct  B-spline  filters  for  the  first  ten  orders. 

Using  Equation  (18)  with  an  approximation  Bp+r+i+i(co)  ~  Bp+r(co)Bi(Lu ),  we  can  obtain  a 
filter  bank  implementation  of  the  transform  decomposition.  The  reconstruction  part 
follows  from  Equations  (8),  (9),  and  (11)  through  (13).  Figure  3  shows  a  filter  bank 
implementation  of  the  transform.  Noninteger  shifts  at  scale  1  are  rounded  to  the  nearest 
integer.  Tables  2  through  5  list  impulse  responses  of  the  filters  used  in  the  filter  bank  for 
p  G  {0,1,2}. 
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Figure  3:  Filter  bank  implementation  of  multiscale  spline  derivatives  for  m  G  [0,  M  —  1]:  (a) 
Prefiltering,  (b)  postfiltering,  (c)  decomposition  module,  and  (d)  reconstruction  module. 
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Table  2:  Impulse  responses  gl{n )  and  g2(n). 


n 

9l{n) 

9\n) 

-1 

1 

1 

0 

-1 

-2 

1 

1 

Table  3:  Impulse  res 


e  res 

ponses 

h(n ),  fcx(n),  k 2 

'■ n ),  ar 

n 

h(n) 

A:1(n) 

t(n) 

-1 

0 

m 

-0.25 

-0.25 

EH 

1 

m 

0.25 

0.25 

Table  4:  Impulse  responses  h(n),  k1( n ),  fc2(n),  and  t(n)  for  p=l. 


n 

h(n ) 

/c1(n) 

k2(n) 

t(n) 

-2 

-1 

B 

-0.0625 

m 

0 

Hi 

-0.3125 

-0.375 

0.375 

1 

0.25 

0.3125 

-0.0625 

0.25 

2 

0.0625 

0.0625 

Table  5:  Impulse  responses  h(ri),  fcx(n),  k2{n ),  and  t(n)  for  p=2. 


n 

h(n ) 

kl(n) 

k2(n) 

t{n) 

H 

-0.015625 

-0.015625 

0.015625 

0.09375 

hi 

-0.109375 

-0.125 

0.234375 

0 

0.375 

-0.34375 

-0.46875 

0.3125 

1 

0.125 

0.34375 

-0.125 

0.234375 

2 

0.109375 

-0.015625 

0.09375 

3 

0.015625 

0.015625 
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Figure  4:  Spline  derivatives  in  the  rr-axis  direction,  (a)  Wavelet  equal  to  the  first  derivative 
of  a  quartic  spline,  (b)  Wavelet  equal  to  the  second  derivative  of  a  quintic  spline. 

The  derived  transform  enables  both  second  derivative  directional  analysis  and  Laplacian  of 
Gaussian  approximations  across  dyadic  scales  (the  latter  can  be  obtained  through 
summation  of  the  outputs  from  blocks  G2(2moj)  applied  along  x  and  y  axis).  Furthermore, 
addition  of  a  block  (?!_.,  (2mcj)  at  each  level  of  decomposition  allows  first  derivative 
directional  analysis  as  well.  Figure  4  shows  first  and  second  derivative  wavelets  obtained  as 
linear  combinations  of  cubic  B-splines. 

2.1.3  Filter  Implementations 

Since  all  two-dimensional  filters  used  in  the  filter  bank  implementation  of  the  transforms 
are  x-y  separable,  only  one-dimensional  filters  need  to  be  implemented.  We  describe  the 
implementation  of  finite  impulse  response  (FIR)  filters  first  and  then  treat  prefiltering  and 
postfiltering  infinite  impulse  response  (HR)  filters  from  Figure  3(a)  and  (b). 

Let  us  refer  to  filters  applied  at  scale  2m  as  filters  at  level  ra+1,  and  let  filters  at  level  1 
(Equations  (10),  (14),  (15),  and  (16)  be  called  “original  filters,”  to  distinguish  them  from 
their  upsampled  versions.  Let  us  split  the  input  to  the  filter  bank  from  Figure  3  into  image 
matrix  rows  and  columns,  each  corresponding  to  a  real  signal  s(n )  G  l2(Z ),  n  G  [0, N  —  1]. 
Depending  on  the  length  of  each  filter  impulse  response,  filtering  an  input  signal  may  be 
computed  either  by  multiplying  the  discrete  Fourier  transforms  of  the  two  sequences  or  by 
circularly  convolving  s(n )  with  a  filter’s  impulse  response.  Using  circular  rather  than  linear 
convolution,  as  is  customary  in  image  processing,  can  lead  to  boundary  artifacts  caused  by 


abrupt  changes  in  the  periodically  extended  signal.  A  common  remedy  for  such  a  problem 
is  realized  by  constructing  a  mirror  extended  signal  [14] 


f  s(—n  —  1)  if  n  G  [-N,  —1] 

|  s{n)  if  n  G  [0,  N  —  1], 


(19) 


where  we  chose  the  signal  sme(n)  to  be  supported  in  [-N,  N  —  1], 

Let  us  classify  symmetric/antisymmetric  real  even-length  signals  into  four  types  [15]: 


Type  I  f(n)  =  f(-n), 


Type  II  /(n)  =  f(-n  -  1), 


Type  III  f(n)  =  -f(-n), 


Type  IV  f(n)  =  -f(-n  -  1), 

where  n  e  [—N,N  —  1].  Note  that  for  Type  I  signals  the  values  at  /( 0)  and  f(—N )  are 
unique,  and  that  for  Type  III  signals  the  values  at  /( 0)  and  f(—N )  are  equal  to  zero. 

(This  is  important  for  storage  requirements:  for  signals  of  Type  II  or  Type  IV,  N  samples 
need  to  be  saved,  while  Type  I  and  Type  III  signals  require  N+l  and  N— 1  sample 
representations,  respectively.) 

Using  properties  of  the  Fourier  transform,  it  is  easy  to  show  that  the  convolution  of 
symmetric/antisymmetric  real  signals  results  in  a  symmetric/antisymmetric  real  signal.  If 
a  symmetric/antisymmetric  real  signal  has  an  even  length,  then  there  always  exists  an 
integer  shift  such  that  the  shifted  signal  belongs  to  one  of  the  above  types. 

Now,  we  are  ready  to  examine  the  filter  bank  implementation  of  the  wavelet  transform 
from  Figure  3  with  filters  given  by  Equations  (10),  (14),  (15),  and  (16)  driven  by  mirrored 
signals  of  the  form  sme(n)  from  Equation  (19)  at  the  input.  Let  the  number  of  levels  M  be 
restricted  by 

M  <  1  +  log2  j-N  *  ,  (20) 

where  Lmax  is  the  length  of  the  longest  original  FIR  filter  impulse  response. 

Each  FIR  filter  block  in  the  filter  bank  consists  of  a  filter  and  a  circular  shift  operator. 
Equation  (20)  guarantees  that  the  length  of  the  filter  impulse  response  does  not  exceed  the 
length  of  the  signal  at  any  block. 

Since  our  mirror  extended  input  row  or  column  sme(n)  is  of  Type  II  and  noninteger  shifts 
at  level  1  are  rounded  to  the  nearest  integer,  it  follows  that  a  processed  one- dimensional 
signal  at  any  point  in  the  filter  bank  belongs  to  one  of  the  types  defined  above.  This  means 
that  filtering  a  signal  of  length  2N  can  be  reduced  to  filtering  a  signal  of  approximately  one 
half  of  its  length. 
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Implementation  is  particularly  simple  for  FIR  filters  designed  with  d= 2  and  p  odd.  Filters 
are  of  Type  I  in  this  case,  so  their  output  will  be  of  Type  II.  An  FIR  filter  block  from  the 
filter  bank  shown  in  Figure  3  can  therefore  be  implemented  by 

L-l 

Fs,mu{n )  =  /(O)ttjj(ra)  +  ^  f(i)[un(n  -  2 mi)  +  un(n  +  2mi)],  n  6  [0,  JV  -  1],  (21) 

i= 1 

where 

u{— n  —  1)  if  n  G  [—  y ,  —1] 

uu(n)  =  <  u\n)  if  n  G  [0,  N  —  1]  (22) 

.  u{2N-n-  1)  if  n  G  [IV,  2A], 

u(n )  is  an  input  signal  to  a  block,  /(n)  is  an  impulse  response  of  G2(2mco),  K2(2mu)), 
T(2moj),  or  H{2muj)  with  p  odd,  L  is  the  length  of  the  filter,  and  N  is  the  length  of  an 
input  signal  s(n)  to  the  filter  bank.  Implementation  of  filters  bp(n)  used  for  prefiltering  and 
postfiltering  (Figure  3(a)  and  (b))  represents  a  special  case  of  Equation  (21)  with  m= 0. 

A  filter  bank  with  the  above  implementation  of  blocks  and  signal  s(n)  at  the  input  yields 
equivalent  results  as  circular  convolution  of  input  sme(n)  as  defined  by  Equation  (19).  In 
addition  to  requiring  one  half  the  amount  of  memory,  the  computational  savings  over  a 
circular  convolution  implementation  of  blocks  are,  depending  on  the  original  filter  length, 
three  to  four  times  fewer  multiplications  and  one  half  as  many  additions. 

A  similar  approach  is  used  for  other  filters.  The  problem  becomes  slightly  more  involved  in 
this  case,  because  the  filters  change  type  from  first  to  subsequent  levels,  and  the  signal 
component  type  can  be  altered  by  a  filter  block  as  well.  As  a  consequence,  an 
implementation  of  blocks  that  use  distinct  original  filters  may  not  be  the  same,  and  the 
implementation  of  blocks  at  level  1  may  differ  from  the  implementation  of  blocks  at  other 
levels  of  analysis. 

The  decomposition  blocks  at  level  1  can  be  implemented  by 

i-i 
2  1 

G-s,ou(n )  =  ]T  g(i)[uH(n  -  i  -  1)  -  un{n  +  i)],  n  G  [1,  N  -  1] 

i= 0 

and 

2  1 

H-Sfiu(n)  =  h(i)[un(n  -  i  -  1)  +  un(n  +  £)],  n  G  [0,  N], 
i= o 

for  p  even,  where  un{l)  is  defined  by  (22),  g(n)  and  h(n)  are  impulse  responses  of  the 
filters  computed  from  (10)  and  (16),  respectively,  and  L  is  the  length  of  the  corresponding 
impulse  response. 

The  output  from  a  block  Gls(cj)  at  level  1  is  of  Type  III,  while  the  output  from  H_s(uj)  at 
the  same  level  is  of  Type  I. 
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The  decomposition  blocks  at  subsequent  levels  m  G  [1,  M— 1]  can  be  implemented  by 

k-i 

G-s,mu(n )  =  £  0(0  -  2m(f  +  0)  -  ui(n  +  2m(*  +  s))]>  n  G  [1,  JV  -  1], 

i=0 

for  p  even, 

k-\ 

G-s,mu(n )  =  9(*)[uii(n  ~  2 m{i  +  a))  -  un(n  +  2 m(i  +  a))],  n  G  [0,  N  -  1], 

i=0 

for  p  odd, 

L—  1 

F-Stmu(n)  =  f(0)ui(n)  +  ^  /(*)[«/(n  -  2mz)  +  u7(n  +  2m«)],  n  G  [0,  TV],  (23) 

t=i 

with  f(n)  =  g(n)  for  d= 2  and  p  even, 

H-SiTnu(n )  =  h(i)[tt/(n  -  2m(z  +  a))  +  u/(n  +  2m(z  +  a))],  n  G  [0,  TV],  (24) 

i=0 

for  p  even,  where 

u(-n)  if  ne  [-y,  -1] 

ui(n)  —  <  u(n)  if  n  G  [0,  N]  (25) 

^u{2N-n)  ifnG  [V  +  l,^]. 

The  outputs  from  blocks  G_s( 2mui)  are  of  Type  III  for  d=l  and  p  even,  of  Type  IV  for 
d=  1  and  p  odd,  and  of  Type  I  for  d= 2  and  p  even,  whereas  the  outputs  from  H-S(2muj) 
are  of  Type  I  for  p  even. 

Next,  the  reconstruction  blocks  at  level  1  can  be  implemented  by 

L 

2 

Klfiu(n)  =  Yj  Hi) [um(n  -i  +  1)  -  uHI(n  +  *)],  n  G  [0,  N  -  1] 

i=l 

and 

L 

2 

HSflu(n)  =  Y  h(i)[ui(n  -  i  +  1)  +  u/(n  +  z)],  n  G  [0,  N  -  1], 

i= 1 

for  p  even,  where 

(  -u(-n)  if  n  G  [— y,  — 1] 

0  if  n  =  0 

um(n )  =  i  u(n)  if  n  G  [1,  N  -  1]  (26) 

0  if  n  =  N 

—u(2N  —  n)  if  n  G  [N  +  1,  ^], 

Ui(n)  is  as  defined  by  (25)  and  k(n)  is  an  impulse  response  of  the  filter  from  (14).  Note 
that  both  outputs  from  blocks  K]{uj)  and  Hs(u)  are  of  Type  II. 
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The  reconstruction  blocks  at  subsequent  levels  can  be  implemented  by 


^-1 

K^mu{n)  —  +  l)[u///(n  —  2m(f  +  s ))  —  um(n  +  2 m(i  +  s))],  n  €.  [0, N], 

i=0 

for  d=  1  and  p  even,  (23)  with  f(n)  —  k(n )  for  d—  2  and  p  even, 

Ks,mu(n )  =  ^2  +  1  )[uiv(n  -  2m(i  +  s))  -  u/v(n  +  2m(i  +  s))],  ne[0,JV-  1], 

i= o 

for  d=l  and  p  odd, 

Hs,mu(n)  —  H—Sjjnu(^Ti) , 

for  p  even,  where  um(l)  is  given  by  (26), 

—u(—n  —  1)  if  n  G  [—  y,  —1] 
uiv(n)  =  <  u(n)  if  n  G  [0,  AT  —  1] 

k-u(2JV-n-l)  ifne[JV,^], 


and  H-a>mu(n)  is  specified  by  Equation  (24).  We  observe  that  the  outputs  from  blocks 
Kg(2mu})  and  Hs(2muj),  m  e  [1,  M  —  1],  for  p  even  are  of  Type  I. 

When  we  compare  the  above  implementation  of  blocks  to  circular  convolution  driven  by  a 
mirrored  signal  sme(n )  at  the  input,  we  observe  that  approximately  twofold  less  memory 
space,  three  to  four  times  fewer  multiplications  and  one  half  as  many  additions  are 
required.  (For  Type  I  signals  an  additional  sample  has  to  be  stored  because  two  values  are 
without  a  pair). 

The  implementation  presented  in  this  section  performs  all  operations  in  the  spatial  domain; 
however,  one  could  also  implement  the  structures  shown  in  Figure  3  with  an  input  signal 
sme(n)  in  the  frequency  domain.  For  short  filter  impulse  responses,  such  as  those  given  in 
Tables  3,  4  and  5,  the  spatial  implementation  described  in  this  section  is  certainly  more 
efficient.  For  long  filter  impulse  responses,  however,  filtering  is  faster  if  implemented  in  the 
frequency  domain.  Additional  details  on  alternative  FIR  filter  implementation  strategies 
can  be  found  in  [16]. 

Implementation  of  HR  filters  b~l{n)  used  for  prefiltering  and  postfiltering  is  a  bit  more 
involved  than  the  one  of  their  FIR  counterparts.  Fortunately,  the  number  of  different  cases 
is  much  smaller  here:  possible  input  to  &pX(n)  in  the  filter  bank  from  Figure  3  is  either  of 
Type  II  or  of  Type  I  (symmetry  types  for  HR  filters  slightly  differ  from  those  defined  for 
FIR  filters:  here,  mirror  extended  signals  are  periodically  repeated,  so  that  they  stretch 
from  —  oo  to  oo).  We  use  ideas  and  a  few  results  from  [17]. 
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Let  us  first  take  a  closer  look  at  the  system  function  Bp  1(z)  with  p  €  {2, 3}.  This  function 
can  be  written  as  a  cascade  of  terms 


E(z)  z_id^l  +  z-i  (l-az-'Xl-az)’ 
which  can  be  expressed  in  a  parallel  form  as 

where  a  and  ^  are  poles  of  the  causal  and  the  anticausal  filter,  respectively. 
The  impulse  response  of  this  term  can  be  written  as 


(27) 


(28) 


We  choose  to  implement  E(z )  in  a  cascade  form  and  therefore  extract  the  difference 
equations  from  Equation  (27): 


c+(n )  =  u(n )  +  ac+(n  —  1)  n  =  1, 2, ,  N—l, 


(29) 


and 

c(n)  =  a  (c(n  +  1)  ~  c+ (n))  n  =  N— 2,  N— 3, . . . ,  0,  (30) 

where  u{ri)  denotes  the  input  to  the  block,  c+(n)  is  the  output  from  the  causal  part,  and 
c(n)  stands  for  the  output  from  the  block. 

To  solve  Equations  (29)  and  (30)  we  need  boundary  conditions  c+(0)  and  c(N— 1).  We 
derive 


c+(0) 


0  N- 1  i+ 1  ,  2  N-i 

S  a  *“«p(0  =  «(°)  +  S  — r — 

i=-oo  i=0  1  “  a 


io 

~  u(0)  +  53  al+1u(z), 


i=0 


and,  using  parallel  form  (28) 


where 


—a 


N-1  aN-i 


_  fyN+l+i 

c(iV-1)  =  r^(c+(iV_1)  +  S  i_a2v — «(*))  - 


-a 


N-1 


1  —  a2 


(c+(iV— 1)  +  ^  aJV  M*))> 

i=N—  1— io 


uup(n) 


un(n  mod  (2N))  if  n  >  0 

un{— (n  +  1)  mod  (2N))  if  n  <  0, 


un(n)  = 


u{n)  if  n  e  [0,  N  —  1] 

u(2N-n-l)  if  n  G  [N,2N  -  1], 


(31) 


(32) 
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N  is  the  length  of  an  input  signal  to  the  filter  bank,  and  io  <  N—l  is  selected  such  that  a 10 
falls  below  a  predefined  precision  threshold. 

For  orders  p  greater  than  three,  we  implement  B~1(z)  as  a  cascade  of  terms  E(z)  with 
different  a’s.  Note  that  the  output  from  block  E(z)  is  always  of  the  same  type  as  the  input 
to  it. 

2.2  Image  Fusion 

Image  fusion  combines  particular  aspects  of  information  from  the  same  imaging  modality  or 
from  distinct  imaging  modalities  and  can  be  used  to  improve  the  reliability  of  a  particular 
computational  vision  task  or  to  provide  a  human  observer  with  a  deeper  insight  about  the 
nature  of  observed  data.  Whether  it  is  combining  different  sensors  or  extending  the 
dynamic  range  of  a  single  sensor,  the  goal  is  to  achieve  more  accurate  inferences  that  can  be 
achieved  by  a  single  sensor  or  a  single  sensor  setting.  By  fusing  together  processed  sections 
of  images,  a  combined  image  which  is  superior  to  the  sum  of  its  parts  can  be  constructed. 
The  simplest  method  of  fusing  images  is  accomplished  by  computing  their  average.  Such  a 
technique  does  combine  features  from  input  images  in  the  fused  image,  however,  the 
contrast  of  the  original  features  can  be  significantly  reduced.  Among  more  sophisticated 
methods,  multiscale  and  multiresolution  analyses  have  become  particularly  popular. 
Different  pyramids  [18,  19]  and  wavelet-based  techniques  [20,  21,  22,  23]  have  been  applied 
to  this  problem. 

In  Section  2.2.1,  we  compared  an  image  fusion  method  based  upon  the  steerable  dyadic 
wavelet  transform  with  recently  published  fusion  methods  based  upon  the  gradient 
pyramid,  the  orthogonal  wavelet  transform,  and  the  biorthogonal  wavelet  transform. 
Section  2.2.2  then  presents  an  application  of  the  fusion  mechanism  to  the  problem  of 
mammographic  feature  contrast  enhancement. 

2.2.1  Comparison  of  Transforms 
Gradient  Pyramid 

Gaussian  pyramid  [4]  was  used  for  construction  of  a  gradient  pyramid  used  in  [18].  Let  the 
generating  filter  kernel  for  the  Gaussian  pyramid  be 
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where  Wb(nx,ny)  denotes  the  3x3  binomial  filter: 

1  [  i  2  1  ' 

WbfjixjUy)  —  —  2  4  2 

10  [  1  2  1 

Level  m  G  N  of  the  Gaussian  pyramid  for  an  input  image  matrix  s(nx,  ny)  is  then 

Qm^ij^xi'^y)  —  *  Qm—l^{p,xi'^'y)')\.2i 


QqS{tIXi  Hy)  — 

Gradient  pyramid  is  obtained  from  the  Gaussian  pyramid  as 


where 


Wy)  —  dj  ♦  ( Qrn^(j^X)  ^y)  “b  Wy  *  /))> 


di  =  [l  -l], 


H  _  i  o-i 

d2  -  T*  i  o 


d3  —  |  ,  and 


d\  — 


j_  -1  0 

V2  0  1 


An  image  is  reconstructed  from  the  gradient  pyramid  by  converting  the  pyramid  to  the 
Laplacian  and  then  to  the  Gaussian  pyramid.  The  Laplacian  pyramid  is  approximated  as 

)  —  fc'm.siP'xy  ^y)  "b  ^  ^y) ;  (33 


where 


JCm.s(jlxi  Hy)  —  o  'y  '  d\  *  'Urns{jlx ,  Uy). 
8  i= 1 


An  approximation  to  the  Gaussian  pyramid  is  obtained  by 

Qm^ij^xi'^'y)  —  £'mS(P,xjHy)  “b  4w  *  ^Qm-\-l^{Plxt  ^y)')'\2- 

Note  that,  since  the  filters  di  have  the  center  of  symmetry  between  samples,  they  need  to 
be  shifted  for  reconstruction,  and  that,  because  of  the  approximation  (33),  the  gradient 
pyramid  does  not  have  the  perfect  reconstruction  property. 
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Discrete  Wavelet  Transform 

Discrete  wavelet  transform  was  implemented  as  a  perfect  reconstruction  filter  bank  [24].  In 
two  dimensions,  the  transform  is  obtained  by  applying  the  one-dimensional  transform 
separately  along  each  dimension.  Level  m  E  N  of  the  transformed  image  matrix  s(nx,ny ) 
is  therefore 


V\#s(n*,  riy)  =  ((Am-is(nx,  ny)  *  g{nx))i2  *  g{ny)) i2, 
W^s(nx,ny)  =  (( Am-is{nx,ny )  *  g(nx))i2  *  h(ny))i2, 
W^s(nx,ny)  =  ((Am-is(nx,  ny)  *  h(nx))42  *  g(ny))i 2,  and 
Ams(nx,ny)  =  ((Am-is(nx,  ny)  *  h(nx))i2  *  h(ny))l2, 


where 


*A-0S(flx)  Tly')  - 


(34) 


{W^s(nx,ny),W^s(nx,ny),W^s(nx,ny)}me[ltM]  and  AMs(nx,  ny)  are  detail  and 
approximation  coefficients  for  M  levels  of  analysis,  respectively,  and  g(n )  and  h(n )  are  the 
decomposition  filters. 

Reconstruction  of  the  approximation  coefficients  at  the  previous  level  is  given  by 


Am-is{nx,  ny)  =  ((W^s(nx,  ny))^2  *  g{ny)  +  (yV^s(nx,  ny))^2  *  h(ny))t2  *  g(nx)+ 
+((W^s(nI,  ny))v  *  g(ny)  +  (.Ams(nx,  ny))t2  *  h(ny))t 2  *  h(nx), 

with  g(n)  and  h{n)  being  the  reconstruction  filters. 

Two-dimensional  wavelets  associated  with  separable  filter  banks,  such  as  the  one  just 
described,  were  constructed  from  tensor  products  of  two  one-dimensional  multiresolution 
analyses  (wavelets  are  products  of  one-dimensional  wavelets  and  scaling  functions)  [6]. 

Note  also  that,  due  to  oversimplified  initialization  (34),  the  discrete  wavelet  transform  may 
be  a  pretty  poor  approximation  to  samples  of  the  continuous  wavelet  transform. 

We  limited  ourselves  to  FIR  filters  (i.e.,  compactly  supported  wavelets).  In  our 
experiments,  we  used  orthogonal  wavelet  transform  with  DAUB12  wavelet  [6],  and 
biorthogonal  wavelet  transform  with  Bior6.8  wavelet  from  MATLAB  Wavelet  Toolbox. 
Image  fusion  is  performed  in  the  transform  space  by  computing  local  statistics  across  the 
decomposition  scales,  and  reconstructing  from  fused  transform  coefficients.  Typical  size  of 
neighborhood  is  between  a  single  pixel  and  5x5  area  with  some  loss  of  contrast  reported 
for  the  latter  [18].  In  general,  the  chosen  size  of  the  area  represents  a  tradeoff  between 
sharpness  and  introduction  of  artifacts.  In  order  to  achieve  our  goal  of  maximum  contrast 
with  minimum  artifacts,  we  limit  the  neighborhood  to  a  single  pixel  (maximizing  contrast) 
and  try  to  choose  the  most  appropriate  transform  (minimizing  artifacts) . 
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Figure  5:  Phantom  used  for  comparisons  of  different  transforms  for  image  fusion. 


After  the  transform  decompositions  of  images  to  be  fused  is  performed,  the  corresponding 
transform  coefficients  of  the  images  are  combined  according  to  the  fusion  rule  into  a  new 
set  of  transform  coefficients  from  which  the  fused  result  is  reconstructed.  As  a  fusion  rule, 
we  used  the  maximum  magnitude  rule  (at  each  position  and  scale  of  the  transforms  the 
coefficient  with  greatest  magnitude  is  chosen)  for  the  gradient  pyramid,  the  orthogonal 
wavelet  transform  (DAUB12),  and  the  biorthogonal  wavelet  transform  (Bior6.8),  and  the 
maximum  oriented  energy  rule  (at  each  position  and  scale  of  the  transforms  the  coefficient 
with  greatest  local  oriented  energy  is  selected)  for  steerable  wavelet  transform  [21]. 

Our  first  two  experiments  used  the  phantom  shown  in  Figure  5.  Image  matrix  has 
dimensions  512  x  512,  and  fusion  was  performed  between  the  original  and  shifted  phantom. 
First,  the  image  to  be  fused  with  the  one  from  Figure  5  was  generated  by  shifting  the 
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Table  6:  Performance  of  fusion  algorithms  based  on  four  different  transforms. 


Transform 

MSE 

MAE 

Gradient  pyramid 

28.66 

13 

Orthogonal  WT 

0.14 

8 

Biorthogonal  WT 

0.21 

7 

Steerable  WT 

0.05 

6 

phantom  one  sample  vertically  towards  the  top  of  the  image.  The  ideal  result  of  fusion 
should  contain  no  double  lines  or  other  artifacts  (the  distance  between  the  corresponding 
points  in  two  images  is  one  pixel,  so  that  the  algorithm  should  merge  shifted  features  into  a 
single  feature).  Figure  6  shows  the  phantom  with  the  surrounding  area  in  the  fused  images. 
The  resulting  images  from  all  four  transforms  were  clipped  (pixel  values  above  the  upper 
limit  of  the  gray  level  range  were  mapped  to  white)  rather  than  scaled.  Please  note  the 
artifacts  present  when  the  orthogonal  and  biorthogonal  wavelet  transforms  were  used.  The 
latter  produced  a  slightly  better  result,  although  artifacts  due  to  aliasing  and  tensor 
product  representation  made  both  fused  images  unacceptable. 

Our  second  experiment  differed  from  the  first  one  only  in  the  fact  that  the  phantom  from 
Figure  5  was  shifted  by  five  samples.  Here,  the  shift  is  large  enough  that  features  from 
both  images  may  be  present  in  the  fused  image.  The  phantom  with  its  surroundings  in  the 
results  of  fusion  using  the  four  transforms  is  demonstrated  in  Figure  7.  Please  note  that 
the  situation  is  similar  to  the  one  in  the  first  experiment.  Both  orthogonal  and 
biorthogonal  wavelet  transforms  exhibited  obvious  artifacts,  while  the  redundant  gradient 
pyramid  and  steerable  dyadic  wavelet  transform  performed  well. 

The  third  experiment  was  geared  towards  a  quantitative  comparison  of  the  four  transforms 
for  image  fusion.  Similar  to  [22],  we  blur  different  parts  of  the  image  and  then  fuse  them  in 
such  a  way  that  each  blurred  part  is  fused  with  its  original  counterpart.  The  ideal  result  of 
fusion  would  be  the  original  image.  Figure  8  shows  the  original  512  x  512  mammogram, 

100  x  100  area  of  interest,  and  two  blurred  images  to  be  fused.  Mean-square  error  (MSE) 
and  maximum  absolute  error  (MAE)  between  the  result  of  fusion  and  the  original  image 
were  computed  for  all  four  methods.  Table  6  summarizes  the  results.  Let  us  remark  that 
these  results  are  rather  typical;  on  a  variety  of  images,  wavelet  based  methods  were  very 
close,  while  consistently  outperforming  the  gradient  pyramid  according  to  the  two  criteria. 
No  significant  artifacts  were  noticed  in  the  “blurring  experiments.”  By  comparing  the 
extracted  regions  with  the  original  ones,  we  subjectively  rated  the  transforms  used  for 
fusion  as:  (1)  steerable  dyadic  wavelet  transform,  (2)  biorthogonal  wavelet  transform,  (3) 
orthogonal  wavelet  transform,  and  (4)  gradient  pyramid.  Again,  the  differences  between 
different  types  of  wavelet  transforms  were  minor,  whereas  the  gradient  pyramid  lagged 
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Figure  6:  Image  fusion  of  phantoms  shifted  by  one  sample,  (a)  Gradient  pyramid,  (b) 
Orthogonal  wavelet  transform,  (c)  Biorthogonal  wavelet  transform,  (d)  Steerable  dyadic 
wavelet  transform. 


Figure  7:  Image  fusion  of  phantoms  shifted  by  five  samples,  (a)  Gradient  pyramid,  (b) 
Orthogonal  wavelet  transform,  (c)  Biorthogonal  wavelet  transform,  (d)  Steerable  dyadic 
wavelet  transform. 


(a)  (b) 


(c)  (d) 


Figure  8:  (a)  Mammogram  with  area  of  interest  delineated,  (b)  Area  of  interest,  (c)  Image 
from  (a)  with  left  half  blurred,  (d)  Image  from  (a)  with  right  half  blurred. 
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behind.  Regions  of  interest  corresponding  to  the  one  from  Figure  8(b)  are  shown  in  Figure 
9.  Wavelet  based  methods  produced  results  that  are  visually  close  to  the  original  region, 
gradient  pyramid,  however,  caused  loss  of  sharpness. 

2.2.2  Fusion  of  Enhanced  Features 

Existing  methods  of  mammograplnc  image  enhancement  can  be  divided  roughly  into  two 
categories:  (1)  methods  aimed  at  better  visualization  of  all  features  present  in  the  image 
[25,  26,  27,  28],  and  (2)  methods  that  target  specific  features  of  importance  (e.g., 
microcalcifications  [29,  30],  stellate  lesions  [11]). 

Methods  from  the  first  category  are  not  optimized  for  a  specific  type  of  cancer  and 
frequently  not  even  for  mammography.  Rather,  they  try  to  improve  the  perceptual  quality 
of  the  entire  image  and  are  often  developed  with  a  framework  more  general  than 
mammography  alone  in  mind. 

The  second  category  methods  concentrate  on  revelation  of  particular  signs  of  malignancy. 
They  can  be  very  successful  in  their  area  of  specialization;  however,  in  order  to  process 
mammogram  for  presence  of  various  features,  one  would  need  to  apply  different  algorithms 
independently  resulting  in  both  larger  number  of  images  to  be  interpreted  by  a  radiologist 
and  increased  computational  complexity  of  such  a  procedure. 

Here,  we  present  an  approach  which  overcomes  these  shortcomings  and  problematic 
limitations  via  synthesis  of  the  two  paradigms  by  means  of  image  fusion. 

The  goal  of  our  method  is  to  adapt  specific  enhancement  schemes  for  distinct 
mammographic  features,  and  then  combine  the  set  of  processed  images  into  an  enhanced 
image.  The  input  mammographic  image  is  first  processed  for  enhancement  of 
microcalcifications,  masses,  and  stellate  lesions.  From  the  resulting  enhanced  images,  the 
final  enhanced  image  is  synthesized  by  means  of  image  fusion.  Wavelet  based  image 
enhancement  and  fusion  are  merged  into  a  unified  framework,  so  that  there  is  no  need  for 
carrying  out  the  two  operations  independently  (i.e.,  computing  wavelet  decompositions, 
modifying  wavelet  coefficients  for  enhancement  of  specific  features,  reconstructing  the 
enhanced  images,  performing  wavelet  transforms  of  the  enhanced  images,  fusing  transform 
coefficients,  and  obtaining  the  final  result  by  reconstruction  from  fused  wavelet 
coefficients).  Both  enhancement  and  fusion  are  therefore  implicit  (i.e.,  performed  in  the 
wavelet  domain  only).  Figure  10  presents  a  block  scheme  of  the  overall  algorithm. 

The  algorithm  consists  of  two  major  steps:  (1)  wavelet  coefficients  are  modified  distinctly 
for  each  type  of  malignancy;  (2)  the  obtained  multiple  sets  of  wavelet  coefficients  are  fused 
into  a  single  set  from  which  the  reconstruction  is  computed.  The  devised  scheme  allows 
efficient  deployment  of  an  enhancement  strategy  appropriate  for  clinical  screening 
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(c)  (d) 


Figure  9:  Region  of  interest  corresponding  to  the  one  from  Figure  8(b)  extracted  from  the 
fused  images  using:  (a)  gradient  pyramid,  (b)  orthogonal  wavelet  transform,  (c)  biorthogonal 
wavelet  transform,  and  (d)  steerable  dyadic  wavelet  transform. 
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Figure  10:  Overview  of  the  algorithm. 


protocols:  enhancement  algorithm  is  first  developed  for  each  specific  type  of  feature 
independently,  and  the  results  are  then  combined  using  an  appropriate  fusion  strategy. 

The  structure  of  the  algorithm  also  enables  independent  development  and  optimization  of 
enhancement  strategies  for  individual  mammographic  features  as  well  as  the  fusion  module. 


Microcalcifications 

Microcalcifications  appear  on  mammograms  in  approximately  half  of  breast  cancer  cases. 
The  assessment  of  shape,  number,  and  distribution  of  microcalcifications  is  important  for  a 
radiologist  to  reach  the  correct  diagnosis.  Microcalcifications  are  smaller  than  1  mm  in  size 
and  can  be  difficult  to  locate  when  they  are  superimposed  on  dense  breast  tissue. 

Several  techniques  have  been  developed  to  improve  the  visibility  of  microcalcifications 
[29,  30,  31,  32].  The  approach  devised  by  Strickland  and  Hahn  [29]  is  particularly  well 
suited  for  our  framework:  they  used  an  undecimated  wavelet  transform  to  approximate 
second  derivatives  of  a  Gaussian  probability  density  function  for  a  multiscale  matched 
filtering  for  presence  of  microcalcifications. 

Strickland  and  Hahn  based  their  method  on  the  observation  that  the  average 
microcalcification  can  be  modeled  by  a  circularly  symmetric  Gaussian  function.  Using  a 
combination  of  a  separable  Markov  process  with  autocorrelation  rnn  =  o^e-QVlfcl+KI  and  a 
nonseparable  Markov  process  with  autocorrelation  rnn  =  a2e~as/k2+l2  to  represent 
mammogram  texture,  they  obtained  the  separable  matched  filter 


)  —  -M"  (cUa, )  Ad  (ojy  ) , 
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and  the  nonseparable  matched  filter 
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Mnsep(ujx,uy)  ~  —  (w*  +  u$e 


T2,,2 


(36) 


In  order  to  deal  with  different  sizes  of  microcalcifications,  one  must  vary  a  of  matched 
filters  (35)  and  (36)  appropriately.  In  Strickland  and  Hahn’s  scheme,  a  wavelet 
decomposition  was  chosen  to  approximate  the  matched  filters  across  the  desired  scale 
range.  Considering  the  lOO^m  resolution  of  the  Nijmegen  database,  the  wavelet  transform 
was  computed  over  the  first  4  octaves.  For  a  denser  sampling  of  scale,  voices  were  inserted 
at  octaves  “2.5”  and  “3.5.”  The  wavelet  analysis  stage  acted  as  a  bank  of  matched  filters; 
wavelet  coefficients  at  locations  indicating  microcalcifications  were  multiplied  by  a  gain 
factor,  and  then  the  inverse  wavelet  transform  was  applied  to  the  modified  coefficients. 

In  our  approach,  microcalcifications  are  modeled  by  central  B-splines.  Using  the  relation 
between  the  standard  deviation  of  a  Gaussian  function  and  the  order  of  B-splines  that 
approximate  it  a  =  \J^ ~  [33],  the  assumption  that  a  Gaussian  object  is  visible 
approximately  over  ±<x  pixels  [29],  and  the  fact  that  the  mammograms  in  the  University  of 
Florida  database  were  digitized  at  116/im  resolution,  four  levels  of  the  transform  from 
Section  2.1  with,  for  example,  p— 3  are  needed  to  encompass  different  sizes  of 
microcalcifications.  The  wavelet  decomposition  including  voices  at  scales  3  and  6 
(corresponding  to  Strickland  and  Hahn’s  octaves  “2.5”  and  “3.5”)  can  be  obtained  by 
deriving  a  counterpart  to  Equation  (18)  for  the  two  scales. 

/3p(3u>)  can  be  related  to  Pp{oj)  by  expressing  /3p(3tu)  as  (cf.  Proposition  1  of  [34]) 


4>(  3w) 


1  ( gin  (tQV+1  / sin(f)\ 

3P+1\sin(| ))  {  |  ) 


Using  Z)m=oe^mw+0i  =  Sm^.  7^t  ^  we  get 

Slnt  2  / 


Pp{  3w)  =  V{u)Pp(u}), 


where 

/I  \  p+1 

V(u,)  =  (e-»  +  1  +  e*)J  . 

Filter  V(co)  can  be  implemented  as  a  moving  sum  with  2(p+l)  additions  per  sample  and  a 
multiplicative  factor  applied  to  the  wavelet  coefficients  [34]. 

Next,  4p(6w)  is  expressed  by  means  of  Equation  (17): 

/3P(6w)  =  H-s(  3u)l3p(3u) 


36 


with  odd  orders  p  used  for  the  sampling  shift  s  to  be  zero. 
Now,  we  can  write 


3~{  ^3S('C)  v)  }  —  S{u!x,uiy)  Br  (u)x)  Br  (u>y)  Bp+r+i+i(u>x)- 

X—nx  ,y=Tly 

■iW KMtiW  Gt,‘( 3u,) 


(37) 


and 


J-{Wls{x,y)  _  _  }  —  S(ux, ujy)  Br  (ujx)  Br  (toy)  Bp+r+i+i(ujx)‘ 

X — Tlx  j2/ — 

■Bp+r+i_H1(w»)Gl7(6a.I)G-:s(6W!,)J7^(3a)l)J/£t',-i(3^1,)V>’+i(1c.«)^-i(«#) 


(38) 


with  notation  being  the  same  as  in  Equation  (18),  and  superscript  p  in  H^s(u>)  denoting 
the  order  p  in  Equation  (16). 

Wavelet  coefficients  obtained  via  Equations  (37)  and  (38)  are  not  used  for 
reconstruction — the  inverse  transform  is  carried  out  as  given  in  Section  2.1. 

The  decomposition  described  by  Equations  (18),  (37),  and  (38)  with  additional  filtering  by 
G2(Ilj)  at  each  scale  l  e  {1, 2, 3, 4, 6, 8}  enables  approximations  to  the  second  derivatives  of 
Gaussian  along  both  x  and  y  directions  and  to  Laplacian  of  Gaussian  across  distinct  scales 
employed  by  Strickland  and  Hahn  [29]  (cf.  Equations  (35)  and  (36)).  We  proceed  in  a 
similar  fashion  as  therein:  the  two  outputs  per  scale  are  thresholded  independently,  all 
binary  results  are  then  combined,  a  circular  region  centered  at  detected  pixel  locations  are 
next  multiplied  by  a  gain,  and,  finally,  the  reconstruction  process  uses  modified  transform 
coefficients. 

Circumscribed  Masses 

Almost  half  of  missed  cancers  appear  on  mammograms  as  masses.  Perception  is  a  problem 
particularly  for  patients  with  dense  fibroglandular  patterns.  The  detection  of  masses  can 
be  especially  difficult  because  of  their  small  size  and  subtle  contrast  compared  with  normal 
breast  structures. 

Fan  and  Laine  [26]  developed  a  discrete  dyadic  wavelet  transform  based  algorithm  suitable 
for  enhancement  of  masses.  They  constructed  an  approximation  to  Laplacian  of  Gaussian 
across  dyadic  scales  for  an  isotropic  input  to  a  piecewise  linear  enhancement  function. 
Approximation  to  Laplacian  of  Gaussian  across  dyadic  scales  is  easy  to  obtain  using 
multiscale  spline  derivatives  derived  in  Section  2.1:  Equation  (18)  with  i  =  0  and  i  —  2 
approximates  the  second  derivative  of  a  Gaussian  function  along  directions  of  x  and  y  axis, 
respectively  (the  corresponding  branches  in  Figure  3(c)  are  the  first  and  third  from  the 
top).  The  appropriate  transform  coefficient  at  each  dyadic  scales  are  therefore  added  and 
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Figure  11:  Enhancement  function  (Equation  (39)  with  K  —  20  and  T  =  1). 

their  sum  input  to  the  piecewise  linear  function 

'  x  -  {K  -  1)T  if  x  <  -T 

C(x)  =  Kx  if  |x|  <  T  (39) 

k  x  +  (K  -  1)T  if  x  >  T. 

used  at  each  level  m+1  of  the  transform  separately.  Due  to  the  expected  size  of  masses, 
levels  greater  than  4  are  enhanced  more  aggressively. 

Figure  11  shows  the  enhancement  function  from  Equation  (39)  for  parameter  values 
K  =  20  and  T  =  1. 

The  multiplicative  factor  obtained  as  the  ratio  between  the  output  and  input  of  the 
enhancement  function  is  next  applied  to  the  original  wavelet  coefficients  [26],  and  then  the 
reconstruction  (Figure  3(b)  and  (d))  is  carried  out. 

Figure  12  shows  the  cranio-caudal  view  of  a  patient  with  bloody  nipple  discharge.  On  the 
enhanced  image  cropped  to  the  area  of  interest,  irregular  anterior  borders  of  a  mass  are 
better  seen. 

Stellate  Lesions 

It  is  important  for  radiologists  to  identify  stellate  lesions  since  their  presence  is  a  serious 
indicator  of  malignancy.  Stellate  lesions  vary  in  size  and  subtlety  and,  in  addition,  do  not 
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(a)  W 

Figure  12:  Contrast  enhancement  of  the  cranio-caudal  of  a  patient  with  bloody  nipple  dis¬ 
charge.  (a)  The  original  mammogram  with  area  of  interest  delineated,  (b)  Unprocessed 
extracted  area,  (c)  The  enhanced  area  improves  the  visualization  of  the  mass. 
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have  a  clear  boundary,  making  them  difficult  to  detect. 

In  the  development  of  our  algorithm,  we  follow  an  observation  made  by  Kegelmeyer  et  al. 
about  the  distortion  of  edge  orientation  distribution  induced  by  a  stellate  lesion  [35]. 
Normal  mammograms  show  a  roughly  radial  pattern  with  structure  radiating  from  the 
nipple  to  the  chest  wall.  A  stellate  lesion  not  only  changes  this  pattern,  but  also  creates 
another  center  from  which  rays  radiate.  Directional  analysis  using  the  Sobel  edge  operator 
was  employed  for  assessment  of  local  orientations  [35]. 

Wavelet  transform  from  Section  2.1  enables  directional  analysis  as  well.  By  adding 
additional  filter  G]_s(2muj)  to  each  scale  of  decomposition  from  Figure  3,  approximations  to 
both  first  and  second  steerable  derivatives  of  a  Gaussian  are  available.  A  multiscale 
derivative-pair  quadratic  feature  detector  is  computed  by  finding  the  maximum  of  the  local 
oriented  energy  with  respect  to  angle  9 

Elm  (x,  y)  =  J(Wle2ms(x,y)y  +  {W2e2ms{x,y))\  (40) 

where  Wld2ms(x,y)  and  W2°2ms(x,  y)  denote  wavelet  decompositions  using  first  (Equation 
(6)  with  d— 1)  and  second  (Equation  (6)  with  d=  2)  derivative  wavelet,  respectively, 
steered  to  angle  9.  The  angle  that  maximizes  the  local  oriented  energy  (40)  represents 
orientation  at  pixel  location  (x,  y). 

Similarly  to  the  method  for  enhancement  of  microcalcifications,  processing  is  carried  out 
within  windows  with  scale  dependent  sizes:  1-norm  of  differences  between  the  local  and 
average  orientations  is  computed  in  the  window  and  used  as  a  measure  of  orientation 
nonuniformity.  Soft  thresholding  as  a  function  of  the  orientation  nonuniformity  measure  is 
next  applied  to  the  transform  coefficients  at  each  dyadic  scale  independently  [11].  The 
altered  coefficients  are  then  included  for  reconstruction. 

Figure  13  shows  the  oblique  view  with  a  mass  visible  in  the  mid-posterior  breast.  The 
enhanced  image  demonstrates  the  irregularity  and  spiculation  of  the  mass. 

Image  fusion  of  the  outputs  from  the  three  modules  is  implicit  rather  than  explicit:  wavelet 
transform  coefficients  are  first  modified  for  enhancement  of  microcalcifications, 
circumscribed  masses,  and  stellate  lesions,  and  then  the  new  coefficients  are  obtained  by 
fusion  before  the  reconstruction  is  accomplished. 

Note  also  that  it  is  possible  to  put  different  weights  on  features,  and  exclude  certain 
features  from  the  final  result. 

2.3  Evaluation 

In  order  to  evaluate  the  contrast  enhancement  ability  of  the  developed  scheme,  we  compare 
the  performance  of  the  wavelet-based  algorithm  described  in  previous  sections  to  results  of 
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Figure  13:  Contrast  enhancement  of  the  oblique  view,  (a)  The  original  mammogram  with 
area  of  interest  delineated,  (b)  Unprocessed  extracted  area,  (c)  The  enhanced  area  improves 
the  visibility  of  the  mass. 


histogram  equalization  and  unsharp  masking. 

Histogram  is  defined  as  a  plot  of  gray  level  probabilities  versus  gray  level  values  in  an 
image.  The  technique  of  obtaining  a  uniform  histogram  (i.e.,  uniform  gray  level 
distribution)  by  reassigning  gray  level  values  is  known  as  histogram  equalization  [36].  It  is 
attractive  due  to  its  speed  and  simplicity;  however,  quantization  errors  can  result  in 
artifacts  such  as  loss  of  edges  [26]. 

Unsharp  masking  is  based  upon  smoothing  an  image  with  a  lowpass  filter  and  then 
subtracting  a  fraction  of  the  filtered  result  from  the  original  image  [37].  To  better  preserve 
edges  in  the  enhanced  image,  the  median  filter  is  a  popular  choice  for  the  smoothing  filter. 


2.3.1  Image  Phantoms 


We  conducted  a  set  of  experiments  using  computer  simulated  phantoms  similar  to  the  ones 
generated  by  Xing  et  al.  [38]. 

The  phantoms  consisted  of  a  505  x  505  image  matrix  with  8  bits  per  pixel  (i.e.,  256  levels 
of  gray)  divided  into  25  (5  x  5  grid)  101  x  101  pixel  squares.  Twelve  randomly  distributed 
squares  contained  a  truncated  Gaussian  shaped  signal  with  the  peak  signal  intensity  I 
(15  <  I  <  35)  levels  above  the  background  level  and  the  radius  of  the  signal  set  to  25 
pixels.  Random  Gaussian  noise  with  the  mean  value  of  128  and  the  standard  deviation  a 
(20  <  a  <  80)  was  added  to  the  image  matrix.  Figure  14(a)  demonstrates  a  distribution  of 
signal  squares  within  the  phantom,  and  Figure  14(b)  shows  the  phantom  with  noise. 

By  choosing  different  values  of  the  peak  signal  intensity  and  of  the  noise  variance  we  can 
alter  the  visibility  of  the  Gaussian  signals.  In  the  course  of  experiments  with  different  peak 
signal  intensities  and  noise  variances,  we  compared  the  ability  of  the  three  contrast 
enhancement  methods  to  improve  the  visibility  of  signals  embedded  in  noise.  We  adopted 
quantitative  criteria  for  measuring  the  performance  of  contrast  enhancement  algorithms 
employed  by  Xing  et  al.  [38] — the  enhancement  factor  EF  was  given  as  the  ratio  between 
the  output  and  input  contrast-to-noise  ratio 


EF  = 


CNR0 

CNR;’ 


where  the  contrast-to-noise  ratio  CNR  was  defined  as  the  ratio  between  contrast 


J_  ^  (g„  -  Bm) 

Ms  m=i  Pm 


and  noise 


1 

m=  1 


m 
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(a)  (b) 


Figure  14:  Generation  of  image  phantoms,  (a)  A  distribution  of  truncated  Gaussian  signals, 
(b)  A  phantom  obtained  by  adding  noise  to  the  image  from  (a). 

with  Sm  and  Bm  denoting  the  gross  count  over  circular  signal  area  in  signal  and 
background  squares,  respectively,  Pm  standing  for  the  number  of  pixels  in  the  ra-th  signal 
area,  Ms  being  equal  to  the  number  of  squares  containing  the  signal,  and  aBm  meaning  the 
standard  deviation  in  the  circular  area  over  which  Bm  was  computed. 

Phantom  images  were  generated  using  the  sets  of  values  I  =  {15, 25, 35}  and 
a  =  {20, 30, 40, 50, 60, 70, 80},  and  contrast  C  and  noise  N  were  computed  for  both  the 
phantom  images  and  enhanced  images  obtained  via  the  three  methods.  These  numerical 
results  are  summarized  in  Table  7. 

The  behavior  of  the  three  methods  in  this  experimental  setup  was  quite  different. 
Histogram  equalization  proved  to  be  the  most  sensitive  to  the  level  of  noise.  In  general  it 
resulted  in  the  largest  contrast  improvement  among  the  three;  however,  it  also  increased 
the  level  of  noise  the  most.  The  resulting  contrast  was  inversely  proportional  to  the  level  of 
noise,  and,  interestingly,  the  noise  after  enhancement  was  pretty  much  independent  of  noise 
in  the  input  image.  It  is  worth  mentioning  that  in  low  to  moderate  noise  conditions  the 
numbers  in  Table  7  do  not  convey  observers’  perceptual  experience  faithfully;  histogram 
equalization  made  objects  more  visible  while  noise,  although  increased,  was  still  acceptable 
to  viewers. 

Unsharp  masking,  on  the  other  hand,  was  not  able  to  increase  contrast,  but  it  dealt  well 
with  noise.  Except  for  the  lowest  noise  cases,  it  always  resulted  in  less  noise  in  the  enhanced 
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Table  7:  Contrast  enhancement  of  image  phantoms  by  means  of  histogram  equalization  (HE),  unsharp  masking  (UM),  and 


CD 

* 


xf 

xf 

rH 

CM 

05 

rH 

o 

CM 

CO 

co 

CM 

05 

00 

xF 

00 

CO 

CM 

xF 

o 

xF 

LO 

CM 

CM 

O 

co 

o 

05 

O 

F— 

05 

F— 

CO 

o 

F- 

CO 

xF 

CO 

I> 

co 

xF 

rH 

f- 

CO 

CO 

rH 

o 

CO 

xF 

CM 

O 

xF 

o 

CO 

CO 

rH 

CO 

rH 

rH 

o 

O 

oo 

to 

CO 

o 

oo 

o 

co 

00 

O 

LO 

05 

CM 

LO 

LO 

oo 

00 

rH 

xF 

i> 

05 

rH 

CO 

w 

rH 

CM 

CM 

CM 

CO 

CO 

CO 

rH 

rH 

rH 

CM 

CM 

CM 

CM 

d 

rH 

rH 

rH 

rH 

CM 

CM 

Eh 

00 

rH 

rH 

CO 

05 

oo 

F~ 

CO 

co 

05 

If- 

F- 

CO 

00 

o 

rH 

LO 

CO 

o 

05 

oo 

CM 

05 

CO 

LO 

05 

CO 

CO 

CO 

rH 

rH 

rH 

CO 

CM 

CO 

03 

co 

00 

00 

F— 

LO 

F~ 

CM 

CO 

CM 

CO 

CM 

rH 

LO 

t> 

05 

oo 

xf 

CM 

LO 

xF 

CO 

LO 

rH 

F- 

CO 

CO 

o 

rH 

XF 

CM 

LO 

05 

CM 

XF 

rH 

rH 

i> 

CO 

o 

xF 

LO 

co 

co 

co 

co 

co 

CO 

LO 

CM 

05 

00 

05 

O 

CO 

f-: 

CO 

xf 

CM 

CO 

cd 

cd 

xF 

cd 

00 

00 

o 

o 

d 

CM 

CM 

rH 

rH 

rH 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CO 

CO 

co 

Eh 

XF 

05 

xf 

C5 

00 

F- 

CO 

CO 

LO 

CM 

LO 

CO 

LO 

05 

F- 

CO 

co 

LO 

F— 

LO 

CM 

F- 

LO 

CM 

F^ 

oo 

O 

LO 

05 

xF 

I> 

rH 

xF 

F— 

O 

05 

CO 

rH 

OO 

oo 

t- 

LQ 

CM 

rH 

CO 

co 

05 

05 

o 

CM 

F— 

rH 

i^ 

xF 

F— 

oo 

CO 

rH 

rH 

CM 

LO 

r~ 

o 

CM 

xF 

co 

CM 

i> 

t> 

t> 

F— 

CM 

00 

xF 

CM 

F- 

LO 

F- 

00 

05 

rH 

F- 

00 

co 

rH 

F- 

CO 

© 

05 

05 

o 

CO 

CM 

05 

l> 

cd 

LO 

LO 

CO 

CO 

CO 

CO 

CM 

rH 

d 

CM 

rH 

rH 

rH 

rH 

CM 

CM 

rH 

rH 

rH 

rH 

r-H 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

05 

rH 

CO 

CO 

xF 

05 

F- 

rH 

If- 

LO 

00 

CM 

CO 

rH 

rH 

05 

o 

t> 

oo 

oo 

rH 

CM 

CO 

o 

CO 

F— 

LO 

05 

rH 

rH 

05 

CO 

xF 

CM 

05 

o 

05 

00 

rH 

CM 

o 

f- 

f- 

F- 

f- 

CO 

LO 

xF 

CO 

F— 

If- 

CO 

CO 

LO 

xF 

CO 

t> 

CO 

co 

co 

LO 

xF 

f- 

f- 

f- 

f~ 

F- 

i> 

F- 

F- 

F- 

F- 

if- 

t> 

F- 

I> 

F- 

F— 

F— 

F- 

F- 

F- 

F- 

w 

o 

o 

d 

o 

o 

o 

O 

d 

d 

o 

d 

d 

d 

O 

O 

O 

O 

O 

O 

O 

d 

o 

05 

XF 

05 

XF 

oo 

CO 

o 

00 

03 

o 

05 

co 

rH 

o 

05 

CO 

oo 

xF 

rH 

05 

xf 

CM 

rH 

© 

o 

LO 

CO 

xF 

CM 

rH 

05 

XF 

04 

co 

xF 

rH 

CM 

CO 

xF 

CD 

xF 

o 

XF 

o 

o 

05 

LO 

rH 

O 

xF 

oo 

1^ 

LO 

05 

LO 

O 

xF 

xF 

xF 

CO 

CD 

05 

CO 

05 

rH 

rH 

CM 

CO 

o 

CO 

05 

oo 

IF- 

05 

rH 

LO 

co 

05 

CO 

xF 

CO 

F- 

O 

XF 

F- 

xf 

CM 

05 

LO 

rH 

xf 

F- 

CO 

rH 

00 

LO 

o 

xF 

F- 

CO 

rH 

00 

xF 

d 

CM 

CM 

CO 

xF 

xF 

LO 

co 

CM 

CM 

CO 

xF 

XF 

LO 

CO 

CM 

CM 

CO 

xF 

xF 

LO 

co 

05 

00 

OO 

O 

CO 

CM 

F- 

CO 

F'- 

05 

rH 

rH 

oo 

XF 

05 

xF 

o 

rH 

xF 

CO 

xF 

05 

on 

05 

CM 

oo 

rH 

F~ 

XF 

05 

LO 

co 

00 

CM 

F- 

o 

CO 

LO 

co 

CM 

F- 

H 

05 

f- 

CO 

O 

co 

CM 

05 

co 

CO 

rH 

05 

05 

CO 

CO 

F— 

05 

05 

co 

CM 

xF 

xF 

o 

XF 

CM 

co 

LO 

CM 

05 

XF 

oo 

rH 

o 

IF- 

CO 

00 

rH 

rH 

05 

CM 

05 

xF 

CO 

F- 

05 

F- 

co 

co 

co 

LO 

LO 

LO 

CM 

rH 

o 

o 

05 

03 

CM 

CO 

LO 

xF 

xF 

co 

CM 

rH 

rH 

rH 

rH 

rH 

CM 

rH 

rH 

rH 

rH 

rH 

rH 

K) 

00 

xf 

CO 

05 

F— 

CM 

o 

o 

CO 

co 

CO 

CM 

rH 

CM 

co 

co 

xF 

CO 

o 

rH 

LO 

V 

rH 

a) 

rH 

o 

co 

O 

rH 

05 

CO 

F- 

CM 

05 

CO 

CO 

CM 

LO 

CO 

on 

oo 

oo 

05 

Ph 

LO 

LO 

CO 

CO 

CO 

F- 

oo 

rH 

xF 

LO 

co 

CO 

F~ 

a> 

F- 

CM 

xF 

LO 

CO 

F- 

00 

05 

05 

05 

05 

05 

05 

05 

05 

05 

05 

03 

05 

05 

05 

00 

05 

05 

05 

05 

05 

05 

H 

O 

o 

o 

o 

o 

o 

o 

d 

o 

d 

o 

o 

d 

d 

o 

O 

O 

o 

d 

o 

d 

fc) 

CO 

f- 

co 

05 

F- 

CO 

oo 

xF 

F- 

00 

co 

rH 

rH 

o 

oo 

CO 

CO 

o 

co 

F- 

05 

CM 

CM 

co 

F- 

rH 

LO 

o 

rH 

LO 

rH 

CO 

CM 

CO 

rH 

F— 

LO 

a) 

co 

xF 

rH 

1- 

XF 

rH 

05 

XF 

o 

XF 

CO 

CM 

LO 

CO 

CO 

LO 

05 

xF 

F- 

o 

xF 

CO 

oo 

LO 

co 

LO 

rH 

CM 

CM 

CM 

CM 

CM 

CM 

00 

03 

CM 

rH 

o 

rH 

CM 

CO 

F— 

05 

05 

o 

rH 

co 

xf 

xF 

XF 

xF 

xF 

xF 

CO 

CO 

CO 

xF 

xF 

xF 

xF 

CM 

CO 

CO 

CO 

CO 

xF 

xF 

f- 

i> 

t> 

F- 

(V 

I> 

F- 

F— 

IF- 

F~ 

F- 

F- 

F~ 

F- 

F- 

F- 

F'- 

F— 

F- 

F- 

F~ 

fcl 

xf 

CO 

F- 

CO 

CO 

F- 

CM 

05 

CM 

CM 

CM 

05 

F— 

05 

F- 

xF 

oo 

CO 

F- 

LO 

05 

v 

00 

00 

LO 

CM 

F- 

oo 

rH 

05 

F^ 

O 

rH 

CO 

LO 

rH 

oo 

05 

rH 

LO 

rH 

F- 

xF 

r? 

CO 

00 

o 

CO 

rH 

o 

LO 

00 

xF 

05 

OO 

00 

rH 

o 

CO 

rH 

rH 

OO 

F- 

05 

O 

xf 

05 

05 

CO 

00 

rH 

00 

05 

CO 

00 

05 

03 

o 

03 

oo 

LO 

CM 

CO 

00 

05 

05 

LO 

CO 

t> 

XF 

rH 

o 

00 

CO 

05 

05 

cd 

05 

F~ 

XF 

xF 

CO 

rH 

CO 

F- 

CO 

d 

CO 

CM 

rH 

rH 

rH 

rH 

LO 

CO 

CM 

CM 

rH 

rH 

rH 

F— 

LO 

xF 

CO 

CM 

CM 

CM 

> 

CM 

o 

00 

LO 

xF 

xF 

CM 

F- 

O 

00 

LO 

xF 

xF 

CM 

F- 

O 

oo 

LO 

xF 

xF 

CO 

05 

05 

CM 

CO 

CM 

co 

co 

05 

03 

CM 

CO 

CM 

CO 

CO 

05 

05 

CM 

<X> 

CM 

CO 

co 

XF 

CM 

CO 

o 

rH 

co 

co 

xf 

CM 

CO 

o 

rH 

CO 

CO 

xF 

CM 

CO 

o 

rH 

co 

00 

F- 

CO 

o 

CO 

rH 

LO 

00 

F- 

CO 

o 

CO 

rH 

LO 

00 

I> 

CO 

o 

CO 

rH 

LO 

05 

05 

05 

05 

F- 

LO 

rH 

03 

05 

05 

05 

F- 

LO 

rH 

05 

05 

05 

05 

F- 

LO 

rH 

t-H 

CM 

CO 

XF 

LO 

CO 

i> 

rH 

CM 

co 

XF 

LO 

CO 

F- 

rH 

CM 

co 

xF 

LO 

CO 

F~ 

LO 

05 

rH 

CO 

rH 

LO 

CM 

CM 

05 

o 

CM 

03 

I> 

CM 

o 

CO 

LO 

rH 

o 

LO 

o 

CM 

LO 

o 

CO 

o 

CO 

CM 

00 

o 

Is- 

CM 

CO 

rH 

o 

xF 

xF 

o 

rH 

F- 

xF 

xF 

LD 

CM 

00 

CO 

CO 

XF 

xF 

oo 

CO 

05 

i> 

rH 

00 

05 

CM 

05 

o 

CM 

LO 

F^ 

O 

O 

05 

oo 

LO 

rH 

CO 

IF- 

if- 

CO 

xF 

o 

co 

LO 

LO 

xF 

xF 

o 

xF 

LO 

xF 

d 

d 

05 

05 

05 

05 

00 

cd 

d 

CO 

cd 

CO 

LO 

xF 

co 

CO 

CO 

CO 

CM 

rH 

d 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

b 

o 

o 

o 

O 

O 

O 

o 

o 

o 

o 

o 

© 

o 

o 

o 

O 

O 

O 

O 

O 

O 

CM 

CO 

XF 

LO 

co 

F- 

00 

CM 

CO 

xF 

LO 

CO 

i> 

00 

CM 

CO 

xF 

LO 

CO 

I> 

00 

H-H 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

LO 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CO 

CO 

CO 

CO 

co 

co 

co 

44 


image.  The  combination  of  loss  of  contrast  and  marginal  noise  improvement,  however, 
resulted  in  worse  visibility  than  in  the  original  images.  This  is  no  surprise  considering  the 
nature  of  the  experiment:  the  Gaussian  shaped  objects  are  comprised  primarily  of  lower 
frequencies  while  unsharp  masking  was  designed  to  emphasize  higher  frequencies. 

The  wavelet  transform  based  approach  was  the  only  one  that  resulted  in  enhancement 
factors  greater  than  one.  It  was  consistently  improving  contrast  and  generally  lowering  the 
amount  of  noise.  High  noise  cases  were  the  ones  where  this  method  really  shined  while  at 
low  levels  of  noise,  the  noise  was,  as  in  unsharp  masking,  amplified  in  the  enhanced  image. 
The  method  did  no  good  when  the  signal  was  strong  and  noise  marginal  (I  =  35,  a  =  20), 
and  such  conditions  also  produced  the  worse  enhancement  factor  in  histogram  equalization. 
Figure  15  displays  an  example  of  moderately  visible  objects.  It  can  be  observed  that 
histogram  equalization  improved  visibility  despite  its  enhancement  factor  being  below  one, 
unsharp  masking  resulted  in  less  noise  but  less  contrast  as  well,  and  wavelet  transform 
yielded  the  best  results.  The  same  signal  constellation  was  used  in  Figure  16.  Here,  one 
can  hardly  notice  any  signal  activity  in  the  original  or  the  histogram  equalized  image,  while 
the  wavelet  transform  based  method  brings  out  the  signal  areas  nicely. 

Please  note  that  the  probability  density  functions  for  signals  and  noise  in  these  image 
phantom  are  not  only  known  but  also  very  simple.  None  of  the  tested  methods  is  taking 
advantage  of  this  fact;  however,  it  would  be  easy  to  devise  one  and  achieve  the  best  results 
for  this  particular  set  of  experiments.  Such  limiting  the  scope  of  a  method,  of  course, 
would  be  diverting  from  the  original  goal — enhancement  of  mammographic  features. 

2.3.2  Phantoms  Blended  Into  Mammograms 

To  compare  the  performance  of  histogram  equalization,  unsharp  masking,  and  wavelet 
based  enhancement  on  mammographic  images,  we  constructed  three  mathematical  models 
of  phantoms  similar  in  appearance  to  the  ones  used  by  Chang  and  Laine  [39]  and  Laine  et 
al.  [25].  The  models  including  features  found  in  circumscribed  masses,  spiculate  lesions, 
and  microcalcifications  are  shown  in  Figure  17. 

These  models  of  varying  intensity,  size,  and  rotation  were  blended  into  mammograms  to  be 
processed  by  the  three  enhancement  techniques  under  comparison.  Quantitative  metric  for 
assessing  the  performance  of  the  techniques  was  adopted  after  Laine  et  al.  [25]:  a  contrast 
improvement  index 

qjj _ Cprocessecj 

^Original 

was  computed  between  the  processed  and  original  contrast  in  the  rectangular  region  of 
interest  snugly  fitting  a  phantom. 
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Figure  15:  Image  phantom  contrast  enhancement,  (a)  Image  phantom  (I  =  35,  a  —  40).  (b) 
Histogram  equalization,  (c)  Unsharp  masking,  (d)  Wavelet  transform  based  enhancement. 
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(a) 


(b) 


Figure  16:  Image  phantom  contrast  enhancement,  (a)  Image  phantom  with  the  same  signal 
constellation  as  in  Figure  15  (I  =  15,  a  =  60).  (b)  Histogram  equalization,  (c)  Unsharp 
masking,  (d)  Wavelet  transform  based  enhancement. 
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Figure  17:  Mathematical  phantoms  that  were  blended  into  mammographic  images. 


The  contrast  was  defined  as 


C  = 


f  +  b 


with  /  and  b  denoting  foreground  ( i.e .,  phantom  pixels)  and  background  (i.e.,  pixels  in  the 
rectangular  area  surrounding  the  phantom),  respectively. 

Note  that  the  input  (reference)  mammograms  for  this  set  of  experiments  were  obtained  by 
applying  contrast  stretching  to  digitized  mammograms.  The  reason  for  this  is  that  we 
wanted  to  avoid  a  commonly  used  “trick”  in  presenting  the  results  of  enhancement:  original 
mammograms  that  occupy  only  a  portion  of  the  available  dynamic  range  are  left  as  they 
are  while  the  enhancement  algorithm  output  utilizes  the  entire  dynamic  range.  We  level  the 
playing  field  by  forcing  the  same  dynamic  range  upon  both  original  and  enhanced  images. 
The  right  cranio-caudal  view  of  Figure  18  shows  a  circumscribed  mass  phantom  blended 
near  the  upper  left  corner  of  the  image  of  matrix  size  1553  x  2048  while  the  original  and 
enhanced  region  of  interest  of  size  121  x  121  are  displayed  underneath  it. 

On  the  oblique  medio-lateral  view  presented  on  Figure  19  of  the  same  size,  blended 
microcalcifications  phantom  slightly  up  and  left  from  the  center  of  the  film  is  practically 
invisible.  The  81  x  81  region  of  interest  is  used  to  compare  the  outputs  of  the  three 
contrast  enhancement  methods. 

Finally,  Figure  20  presents  an  example  of  a  spiculate  lesion  phantom  enhancement.  The 
phantom  is  too  thin  to  be  visible  in  the  upper  right  corner  of  the  mammographic  image  of 
size  1553  x  2048,  but  it  does  appear  on  the  extracted  81  x  81  region  of  interest  from  the 
original  and  the  enhanced  images. 

Comparison  of  the  enhancement  methods  by  means  of  both  visual  appearance  and 
quantitative  criteria  demonstrated  that  the  wavelet  based  approach  outperformed 
histogram  equalization  and  unsharp  masking.  The  CII  values  given  in  Tables  8,  9,  and  10 
for  different  mammograms,  phantom  intensities,  sizes,  and  orientations  show  that 


Figure  18:  Contrast  enhancement  of  a  mass  resembling  phantom,  (a)  Original  mammogram 
with  a  phantom  blended  near  the  upper  left  corner,  (b)  Region  of  interest  from  the  original 
mammogram,  (c)  Histogram  equalization,  (d)  Unsharp  masking,  (e)  Wavelet  transform 
based  enhancement. 


(a) 


Figure  19:  Contrast  enhancement  of  a  microcalcifications  resembling  phantom,  (a)  Original 
mammogram  with  a  phantom  blended  slightly  up  and  left  from  the  center  of  the  film,  (b) 
Region  of  interest  from  the  original  mammogram,  (c)  Histogram  equalization,  (d)  Unsharp 
masking,  (e)  Wavelet  transform  based  enhancement. 
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Figure  20:  Contrast  enhancement  of  a  microcalcifications  resembling  phantom,  (a)  Original 
mammogram  with  a  phantom  blended  slightly  up  and  left  from  the  center  of  the  film,  (b) 
Region  of  interest  from  the  original  mammogram,  (c)  Histogram  equalization,  (d)  Unsharp 
masking,  (e)  Wavelet  transform  based  enhancement. 


51 


Table  8:  Contrast  enhancement  of  circumscribed  mass  phantoms  by  means  of  histogram 
equalization  (HE),  unsharp  masking  (UM),  and  wavelet  transform  (WT). 


Case 

CII  he 

CII  UM 

Cllwr 

1 

1.6917 

1.4749 

2.4790 

2 

0.6368 

0.5615 

2.5695 

3 

1.1028 

0.5702 

2.9598 

4 

0.7780 

0.5644 

3.5478 

5 

2.1308 

0.4583 

1.5779 

6 

0.6249 

0.6345 

3.0910 

1.9501 

1.2311 

2.6068 

1.4860 

0.8913 

2.7621 

1.4565 

1.0185 

1.8214 

0.9218 

0.7382 

2.1763 

Table  9:  Contrast  enhancement  of  circumscribed  mass  phantoms  by  means  of  histogram 
equalization  (HE),  unsharp  masking  (UM),  and  wavelet  transform  (WT). 


Case 

CIIhe 

CIW 

CIIyVT 

1 

0.5335 

1.3471 

3.0844 

2 

1.1018 

0.9552 

2.4362 

3 

1.8519 

0.8615 

2.3966 

4 

0.9553 

0.9345 

2.3869 

5 

1.4057 

0.9355 

2.9169 

6 

1.4103 

0.8936 

3.0579 

1.3529 

0.8132 

3.0099 

1.1389 

1.2028 

3.1987 

0.6038 

1.2722 

2.7468 

1.4451 

1.4660 

3.4186 

histogram  equalization  outperformed  wavelet  transform  based  enhancement  in  certain 
situations.  This  is  due  to  the  fact  that  the  phantom  has  basically  no  effect  on  the  gray 
level  distribution,  and  its  gray  levels  sometimes  simply  got  quantized  to  high  enough  a 
value  for  a  significant  contrast  improvement.  More  than  inconsistency  of  histogram 
equalization  is  worrisome  its  tendency  to  introduce  artifacts.  This  was  particularly  obvious 
in  enhancement  of  microcalcifications  phantoms;  as  exhibited  in  Figure  19(c),  the  circles 
were  completely  distorted.  Unsharp  masking  was  both  consistent  and  artifact  free,  but 
unfortunately  resulted  in  least  enhancement.  It  seems  that  the  images  were  too  complex 
and  the  phantoms  too  subtle  for  unsharp  masking  to  make  a  difference. 
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Table  10:  Contrast  enhancement  of  spiculate  lesion  phantoms  by  means  of  histogram  equal¬ 
ization  (HE),  unsharp  masking  (UM),  and  wavelet  transform  (WT). 


Case 

CII  HE 

CIIwt 

1 

2.3008 

1.3125 

3.3071 

2 

1.8462 

1.2026 

2.6721 

3 

0.8381 

1.0196 

1.8813 

4 

1.3795 

0.8318 

2.5028 

5 

1.7095 

1.4289 

2.3046 

6 

2.1897 

1.1934 

1.9822 

1.3028 

0.5417 

2.1509 

2.6979 

1.3784 

2.8600 

2.8537 

0.5936 

2.4966 

3 

0.8998 

0.8216 

2.6449 
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3  Key  Research  Accomplishments 

•  A  redundant  wavelet  transform  that  enables  artifact-free  translation  and  rotation 
invariant  processing  was  constructed  and  efficiently  implemented  as  a  bank  of  digital 
filters. 

•  A  modular  scheme  for  fusion  of  locally  enhanced  mammographic  features  was  devised 
such  that  enhancement  algorithms  targeting  specific  features  can  be  further 
developed  and  optimized  separately  from  the  fusion  mechanism. 

•  The  derived  wavelet  transform  proved  to  be  more  suitable  for  image  fusion  than 
traditional  orthogonal/biorthogonal  wavelet  transforms  and  image  pyramids. 

•  Quantitative  evaluation  gave  the  wavelet  based  enhancement  technique  an  edge  over 
histogram  equalization  and  unsharp  masking  while  at  the  same  time  no  harmful 
artifacts  emanating  from  the  wavelet  method  were  observed. 
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4  Reportable  Outcomes 


Publications 

[1]  I.  Koren,  A.  Laine,  S.  Smith,  E.  NickolofF,  and  F.  Taylor,  “Visualization  of 
mammograms  via  fusion  of  enhanced  features,”  in  K.  Doi,  H.  MacMahon,  M.  L.  Giger,  and 
K.  R.  Hoffmann,  Eds.,  Computer-Aided  Diagnosis  in  Medical  Imaging ,  Amsterdam,  The 
Netherlands:  Elsevier,  1999. 

[2]  I.  Koren,  A.  Laine,  and  F.  Taylor,  “Enhancement  via  fusion  of  mammographic  features,” 
in  Proc.  IEEE  Int.  Conf.  Image  Process .,  Chicago,  II,  Oct.  1998,  vol.  1,  pp.  722-726. 

[3]  I.  Koren  and  A.  Laine,  “A  discrete  dyadic  wavelet  transform  for  multidimensional 
feature  analysis,”  in  M.  Akay  (Editor),  Time- Frequency  and  Wavelets  in  Biomedical  Signal 
Engineering ,  New  York,  NY:  IEEE  Press,  1998,  pp.  425-449. 

[4]  I.  Koren,  A.  Laine,  and  F.  Taylor,  “An  overcomplete  enhancement  of  digital 
mammograms,”  Era  of  Hope,  A  Multidisciplinary  Reporting  of  DoD  Progress,  Washington, 
D.C.,  Oct.-Nov.  1997,  vol.  1,  pp.  107-108. 

[5]  A.  F.  Laine,  I.  Koren,  S.  Schuler,  W.  Huda,  and  B.  G.  Steinbach,  “Contrast 
enhancement  of  mammographic  features  via  multiscale  analysis,”  RSNA  82nd  Scientific 
Assembly  and  Annual  Meeting,  Chicago,  IL,  1996. 

Funded  Research 

National  Cancer  Institute,  National  Institutes  of  Health,  “Athena  Mammographic 
Technology,”  1999  to  2000,  $125,000;  Iztok  Koren,  Principal  Investigator. 

Employment 

The  Athena  Group,  Inc.,  Gainesville,  FL.  Iztok  Koren,  development  of  mammographic 
image  manipulation  and  enhancement  system. 
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5  Conclusions 


In  the  course  of  this  project,  we  developed  a  method  that  improves  the  visibility  of  specific 
mammographic  features  by  local  enhancement  and  fusion  of  the  enhanced  areas.  The 
described  method  incorporates  a  variety  of  properties  of  mammographic  image 
enhancement  methods  tailored  to  specific  signs  of  malignancy  into  a  unified  computational 
framework. 

Our  method  is  based  upon  a  wavelet  transform  which  we  constructed  such  that  its  analysis 
stage  enabled  approximations  to  directional  first  and  second  derivatives  of  a  Gaussian 
function  and  to  Laplacian  of  Gaussian  across  distinct  scales.  Such  a  decomposition  is 
suitable  for  both  anisotropic  and  isotropic  multiscale  analysis,  and,  in  addition,  provides  an 
adaptable  framework  for  incorporation  of  a  variety  of  mammographic  processing  methods. 
The  derived  steerable  dyadic  wavelet  transform  has  also  proved  flexible  enough  for 
enhancement  and  fusion  of  individual  types  of  mammographic  features.  Separate 
enhancement  algorithms  have  been  developed  for  microcalcifications,  circumscribed  masses, 
and  stellate  lesions,  and  fusion  of  the  modified  transform  coefficients  performed  before  the 
reconstruction  of  the  final  enhanced  image.  The  devised  algorithm  is  also  well  suited  for 
further  refinements;  optimizations  can  be  performed  for  each  type  of  malignancy  alone,  and 
separately  for  the  fusion  module. 

It  is  of  prime  importance  that  neither  the  enhancement  nor  the  fusion  process  introduce 
any  artifacts  that  could  adversely  impact  the  radiologists’  decisions.  We  compared  the 
steerable  dyadic  wavelet  transform  with  the  gradient  pyramid,  orthogonal  wavelet 
transform,  and  biorthogonal  wavelet  transform  for  fusion  of  features  relevant  to 
mammography.  During  our  experiments,  the  steerable  dyadic  wavelet  transform  was 
exhibiting  a  combination  of  best  properties  of  the  gradient  pyramid  and  of  the 
orthogonal/biorthogonal  wavelet  transform.  The  steerable  dyadic  wavelet  transform 
behaved  similarly  as  the  gradient  pyramid  in  a  sense  that  it  did  not  introduce  artifacts 
commonly  present  when  the  orthogonal  and  biorthogonal  wavelet  transforms  were  used.  At 
the  same  time,  the  steerable  dyadic  wavelet  transform  outperformed  the  gradient  pyramid 
by  acting  like  the  orthogonal  and  biorthogonal  wavelet  transforms  in  terms  of 
mathematical  criteria  and  sharpness  in  the  fused  image. 

Our  final  task  was  evaluation  of  the  developed  wavelet  based  enhancement  method.  We 
adopted  quantitative  criteria  to  measure  the  performance  of  our  contrast  enhancement 
technique  on  simulated  phantoms  embedded  in  noise  and  on  mammographic  feature 
phantoms  blended  into  mammograms.  The  wavelet  based  method  outperformed  histogram 
equalization  and  unsharp  masking  with  respect  to  both  quantitative  criteria  and  absence  of 
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undesirable  artifacts. 

This  work  represents  an  important  step  towards  development  of  a  system  which  can  be 
deployed  in  a  clinical  setting.  In  order  to  get  closer  to  that  goal,  a  receiver  operating 
characteristic  (ROC)  study  which  was  unfortunately  beyond  the  scope  of  this  project  will 
have  to  be  carried  out. 


57 


References 


[1]  S.  H.  Landis,  T.  Murray,  S.  Bolden,  and  P.  A.  Wingo,  “Cancer  statistics,  1999,”  CA 
Cancer  J.  Clin.,  vol.  49,  no.  1,  pp.  8-31,  1999. 

[2]  R.  A.  Smith,  “Epidemiology  of  breast  cancer,”  in  A  Categorical  Course  in  Physics, 
Technical  Aspects  of  Breast  Imaging,  A.  G.  Haus  and  M.  J.  Yaffe,  Eds.  The  79th 
Scientific  Assembly  and  Annual  Meeting  of  the  Radiological  Society  of  North  America 
(RSNA),  1993,  pp.  21-33. 

[3]  M.  J.  T.  Smith  and  T.  P.  Barnwell,  III,  “Exact  reconstruction  techniques  for 
tree-structured  subband  coders,”  IEEE  Trans.  Acoust.  Speech  Signal  Process.,  vol.  34, 
no.  3,  pp.  434-441,  1986. 

[4]  P.  J.  Burt  and  E.  H.  Adelson,  “The  Laplacian  pyramid  as  a  compact  image  code,” 
IEEE  Trans.  Commun.,  vol.  31,  no.  4,  pp.  532-540,  1983. 

[5]  A.  Witkin,  “Scale  space  filtering,”  in  Proc.  Int.  Joint  Conf.  Artif.  Intel! ,  Karlsruhe, 
Germany,  1983,  pp.  1019-1022. 

[6]  I.  Daubechies,  Ten  Lectures  on  Wavelets,  SIAM,  Philadelphia,  PA,  1992. 

[7]  M.  Holschneider,  R.  Kronland-Martinet,  J.  Morlet,  and  Ph.  Tchamitchian,  “A 
real-time  algorithm  for  signal  analysis  with  the  help  of  the  wavelet  transform,”  in 
Wavelets,  Time- Frequency  Methods  and  Phase  Space,  J.  M.  Combes,  A.  Grossmann, 
and  Ph.  Tchamitchian,  Eds.,  Springer- Verlag,  Berlin,  Germany,  1989,  pp.  286-297. 

[8]  J.  Canny,  “A  computational  approach  to  edge  detection,”  IEEE  Trans.  Pattern  Anal. 
Mach.  Intel!,  vol.  8,  no.  6,  pp.  679-698,  1986. 

[9]  W.  T.  Freeman  and  E.  H.  Adelson,  “The  design  and  use  of  steerable  filters,”  IEEE 
Trans.  Pattern  Anal.  Mach.  Intel!,  vol.  13,  no.  9,  pp.  891-906,  1991. 

[10]  J.  Babaud,  A.  P.  Witkin,  M.  Baudin,  and  R.  O.  Duda,  “Uniqueness  of  the  Gaussian 
kernel  for  scale-space  filtering,”  IEEE  Trans.  Pattern  Anal.  Mach.  Intel!,  vol.  8,  no.  1, 
pp.  26-33,  1986. 

[11]  I.  Koren,  A.  Laine,  F.  Taylor,  and  M.  Lewis,  “Interactive  wavelet  processing  and 
techniques  applied  to  digital  mammography,”  in  Proc.  IEEE  Int.  Conf.  Acoust.  Speech 
Signal  Process.,  Atlanta,  GA,  1996,  vol.  3,  pp.  1415-1418. 


58 


t  t 


[12]  I.  Koren,  A  multiscale  spline  derivative-based,  transform  for  image  fusion  and 
enhancement ,  PhD  thesis,  Department  of  Electrical  and  Computer  Engineering, 
University  of  Florida,  Gainesville,  FL,  1996. 

[13]  A.  J.  Jerri,  “The  Shannon  sampling  theorem — its  various  extensions  and  applications: 
A  tutorial  review,”  Proc.  IEEE,  vol.  65,  no.  11,  pp.  1565-1596,  1977. 

[14]  I.  Koren  and  A.  Laine,  “A  discrete  dyadic  wavelet  transform  for  multidimensional 
feature  analysis,”  in  Time- Frequency  and  Wavelets  in  Biomedical  Signal  Engineering, 
M.  Akay,  Ed.,  IEEE  Press,  New  York,  NY,  1997,  pp.  425-449. 

[15]  A.  V.  Oppenheim  and  R.  W.  Schafer,  Discrete-Time  Signal  Processing ,  Prentice-Hall, 
Englewood  Cliffs,  NJ,  1989. 

[16]  O.  Rioul  and  P.  Duhamel,  “Fast  algorithms  for  discrete  and  continuous  wavelet 
transforms,”  IEEE  Trans.  Inf.  Theory,  vol.  38,  no.  2,  pp.  569-586,  1992. 

[17]  M.  Unser,  A.  Aldroubi,  and  M.  Eden,  “Fast  B-spline  transforms  for  continuous  image 
representation  and  interpolation,”  IEEE  Trans.  Pattern  Anal.  Mach.  Intel l,  vol.  13, 
no.  3,  pp.  277-285,  1991. 

[18]  P.  J.  Burt  and  R.  J.  Kolczynski,  “Enhanced  image  capture  through  fusion,”  in  Proc. 
Int.  Conf.  Comput.  Vision ,  Berlin,  Germany,  1993,  pp.  173-182. 

[19]  A.  Toet,  “Image  fusion  by  a  ratio  of  low-pass  pyramid,”  Pattern  Recognit.  Lett.,  vol. 

9,  no.  4,  pp.  245-253,  1989. 

[20]  L.  J.  Chipman,  T.  M.  Orr,  and  L.  N.  Graham,  “Wavelets  and  image  fusion,”  in  Proc. 
IEEE  Int.  Conf.  Image  Process.,  Washington,  D.C.,  1995,  vol.  3,  pp.  248-251. 

[21]  I.  Koren,  A.  Laine,  and  F.  Taylor,  “Image  fusion  using  steerable  dyadic  wavelet 
transform,”  in  Proc.  IEEE  Int.  Conf.  Image  Process.,  Washington,  D.C.,  1995,  vol.  3, 
pp.  232-235. 

[22]  H.  Li,  B.  S.  Manjunath,  and  S.  K.  Mitra,  “Multi-sensor  image  fusion  using  the  wavelet 
transform,”  in  Proc.  IEEE  Int.  Conf.  Image  Process.,  Austin,  TX,  1994,  vol.  1,  pp. 
51-55. 

[23]  T.  Ranchin,  L.  Wald,  and  M.  Mangolini,  “Efficient  data  fusion  using  wavelet 
transform:  the  case  of  SPOT  satellite  images,”  in  Mathematical  Imaging:  Wavelet 
Applications  in  Signal  and  Image  Processing,  Proc.  SPIE,  A.  F.  Laine,  Ed.,  San  Diego, 
CA,  1993,  vol.  2034,  pp.  171-178. 


59 


[24]  S.  G.  Mallat,  “A  theory  for  multiresolution  signal  decomposition:  the  wavelet 

representation,”  IEEE  Trans.  Pattern  Anal.  Mach.  Intel l.,  vol.  11,  no.  7,  pp.  674-693, 
1989. 


*  *  * 


[25]  A.  F.  Laine,  S.  Schuler,  J.  Fan,  and  W.  Huda,  “Mammographic  feature  enhancement 
by  multiscale  analysis,”  IEEE  Trans.  Med.  Imaging,  vol.  13,  no.  4,  pp.  725-740,  1994. 

[26]  J.  Fan  and  A.  Laine,  “Multiscale  contrast  enhancement  and  denoising  in  digital 
radiographs,”  in  Wavelets  in  Medicine  and  Biology,  A.  Aldroubi  and  M.  Unser,  Eds., 
CRC  Press,  Boca  Raton,  FL,  1996,  pp.  163-189. 

[27]  C.-M.  Chang  and  A.  Laine,  “Enhancement  of  mammograms  from  oriented 
information,”  in  Proc.  IEEE  Int.  Conf.  Image  Process.,  Santa  Barbara,  CA,  1997, 
vol.  3,  pp.  524-527. 

[28]  W.  M.  Morrow,  R.  B.  Paranjape,  R.  M.  Rangayyan,  and  J.  E.  L.  Desautels, 
“Region-based  contrast  enhancement  of  mammograms,”  IEEE  Trans.  Med.  Imaging, 
vol.  11,  no.  3,  pp.  392-406,  1992. 

[29]  R.  N.  Strickland  and  H.  I.  Hahn,  “Wavelet  transforms  for  detecting  microcalcifications 
in  mammograms,”  IEEE  Trans.  Med.  Imaging,  vol.  15,  no.  2,  pp.  218-229,  1996. 

[30]  H.  Yoshida,  W.  Zhang,  W.  Cai,  K.  Doi,  R.  M.  Nishikawa,  and  M.  L.  Giger, 
“Optimizing  wavelet  transform  based  on  supervised  learning  for  detection  of 
microcalcifications  in  digital  mammograms,”  in  Proc.  IEEE  Int.  Conf.  Image  Process., 
1995,  vol.  3. 

[31]  L.  W.  Estevez  and  N.  D.  Kehtarnavaz,  “Computer  assisted  enhancement  of 
mammograms  for  detection  of  microcalcifications,”  in  Proc.  IEEE  Symp.  Comput. 
Based  Med.  Syst.,  1995,  pp.  16-23. 

[32]  H.  Li,  K.  J.  R.  Liu,  and  S.  C.  B.  Lo,  “Fractal  modeling  of  mammogram  and 
enhancement  of  microcalcifications,”  in  Proc.  IEEE  Nucl.  Sci.  Symp.  Med.  Imaging 
Conf.,  1997,  vol.  3,  pp.  1850-1854. 

[33]  M.  Unser,  A.  Aldroubi,  and  M.  Eden,  “On  the  asymptotic  convergence  of  B-spline 
wavelets  to  Gabor  functions,”  IEEE  Trans.  Inf.  Theory,  vol.  38,  no.  2,  pp.  864-872, 
1992. 

[34]  M.  Unser,  A.  Aldroubi,  and  S.  J.  Schiff,  “Fast  implementation  of  the  continuous 
wavelet  transform  with  integer  scales,”  IEEE  Trans.  Signal  Process.,  vol.  42,  no.  12, 
pp.  3519-3523,  1994. 


60 


[35]  W.  P.  Kegelmeyer,  J.  M.  Pruneda,  P.  D.  Bourland,  A.  Hillis,  M.  V.  Riggs,  and  M.  L. 
Nipper,  “Computer-aided  mammographic  screening  for  spiculated  lesions,”  Radiology , 
vol.  191,  no.  2,  pp.  331-337,  1994. 

[36]  R.  C.  Gonzales  and  R.  C.  Woods,  Digital  Image  Processing ,  Addison- Wesley,  Reading, 
MA,  1992. 

[37]  J.  S.  Lim,  Two-Dimensional  Signal  and  Image  Processing,  Prentice  Hall,  Englewood 
Cliffs,  NJ,  1990. 

[38]  Y.  Xing,  W.  Huda,  A.  Laine,  and  J.  Fan,  “Simulated  phantom  images  for  optimizing 
wavelet  based  image  processing  algorithms  in  mammography,”  in  Mathematical 
Methods  in  Medical  Imaging  III,  Proc.  SPIE,  F.  L.  Bookstein,  J.  S.  Duncan,  N.  Lange, 
and  D.  C.  Wilson,  Eds.,  San  Diego,  CA,  1994,  vol.  2299,  pp.  207-217. 

[39]  C.-M.  Chang  and  A.  Laine,  “Coherence  of  multiscale  features  for  enhancement  of 
digital  mammograms,”  IEEE  Trans.  Inf.  Technol.  Biomed.,  vol.  3,  no.  1,  pp.  32-46, 
1999. 


61 


List  of  Personnel 

Iztok  Koren,  Ph.D.,  Principal  Investigator 


62 


